{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute health expenditures per income level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats \n",
    "from importlib import reload\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import os\n",
    "import sys\n",
    "from statsmodels.formula.api import ols,wls\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model import micro,macro, params\n",
    "from model import distributions as dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulated Gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "file = open('output/eq_pus_ge3.pkl','rb')\n",
    "eq = pickle.load(file)\n",
    "file.close()\n",
    "opt = pd.read_pickle('output/opt_pus_ge3.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wmean(x,var):\n",
    "    xx = x.loc[~x[var].isna(),:]\n",
    "    names = {var: (xx[var] * xx['ps']).sum()/xx['ps'].sum()} \n",
    "    return pd.Series(names, index=[var])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"10\" halign=\"left\">m</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>e</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>h</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.525576</td>\n",
       "      <td>7.810403e-01</td>\n",
       "      <td>1.29577</td>\n",
       "      <td>2.689467</td>\n",
       "      <td>7.29176</td>\n",
       "      <td>13.616625</td>\n",
       "      <td>16.92276</td>\n",
       "      <td>18.823232</td>\n",
       "      <td>2.058494e+01</td>\n",
       "      <td>22.230655</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>5.656805e-18</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.268590e-14</td>\n",
       "      <td>0.019326</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          m                                                                 \\\n",
       "e         0             1        2         3        4          5         6   \n",
       "h                                                                            \n",
       "0  0.525576  7.810403e-01  1.29577  2.689467  7.29176  13.616625  16.92276   \n",
       "1  0.000000  5.656805e-18  0.00000  0.000000  0.00000   0.000000   0.00000   \n",
       "\n",
       "                                       \n",
       "e          7             8          9  \n",
       "h                                      \n",
       "0  18.823232  2.058494e+01  22.230655  \n",
       "1   0.000000  4.268590e-14   0.019326  "
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab_eh = opt.groupby(['h','e']).apply(wmean,var='m').unstack()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   e\n",
       "h  0    0.677376\n",
       "   1    0.719265\n",
       "   2    0.778303\n",
       "   3    0.846530\n",
       "   4    0.910622\n",
       "   5    0.949463\n",
       "   6    0.962733\n",
       "   7    0.965976\n",
       "   8    0.967188\n",
       "   9    0.967949\n",
       "dtype: float64"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab_h = opt.groupby('e').apply(wmean,var='h').unstack()\n",
    "tab_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>m</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>e</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.169564</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.219266</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.287269</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.412752</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.651721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.688138</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.630657</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.640445</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.675423</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.731230</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          m\n",
       "e          \n",
       "0  0.169564\n",
       "1  0.219266\n",
       "2  0.287269\n",
       "3  0.412752\n",
       "4  0.651721\n",
       "5  0.688138\n",
       "6  0.630657\n",
       "7  0.640445\n",
       "8  0.675423\n",
       "9  0.731230"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_sim = opt.groupby('e').apply(wmean,var='m')\n",
    "grad_m_sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>m</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.739822</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.663559</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.559215</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.366673</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.055878</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>-0.032321</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>-0.017302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.036369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.121999</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           m\n",
       "1  -0.739822\n",
       "2  -0.663559\n",
       "3  -0.559215\n",
       "4  -0.366673\n",
       "5   0.000000\n",
       "6   0.055878\n",
       "7  -0.032321\n",
       "8  -0.017302\n",
       "9   0.036369\n",
       "10  0.121999"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_sim['m'] = grad_m_sim['m']/grad_m_sim.loc[4,'m']-1\n",
    "grad_m_sim.index = [x for x in range(1,11)]\n",
    "grad_m_sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>duid</th>\n",
       "      <th>dupersid</th>\n",
       "      <th>famidyr</th>\n",
       "      <th>famszeyr</th>\n",
       "      <th>famrfpyr</th>\n",
       "      <th>region53</th>\n",
       "      <th>age</th>\n",
       "      <th>racex</th>\n",
       "      <th>hispan</th>\n",
       "      <th>marry53x</th>\n",
       "      <th>...</th>\n",
       "      <th>nh_diabe</th>\n",
       "      <th>nh_lunge</th>\n",
       "      <th>nh_cancre</th>\n",
       "      <th>smokev</th>\n",
       "      <th>nh_smoken</th>\n",
       "      <th>mergeMEPSNHIS</th>\n",
       "      <th>smoken</th>\n",
       "      <th>totinc</th>\n",
       "      <th>inscov</th>\n",
       "      <th>_merge</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002014</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>30000.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002014</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>25200.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002021</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002021</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002038</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 125 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      duid  dupersid famidyr  famszeyr  famrfpyr  region53   age  racex  \\\n",
       "0  20002.0  20002014       A       6.0       1.0       4.0  39.0    1.0   \n",
       "1  20002.0  20002014       A       6.0       1.0       4.0  40.0    1.0   \n",
       "2  20002.0  20002021       A       6.0       0.0       4.0  16.0    1.0   \n",
       "3  20002.0  20002021       A       6.0       0.0       4.0  17.0    1.0   \n",
       "4  20002.0  20002038       A       6.0       0.0       4.0  14.0    1.0   \n",
       "\n",
       "   hispan  marry53x  ...  nh_diabe  nh_lunge  nh_cancre  smokev  nh_smoken  \\\n",
       "0     1.0       1.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "1     1.0       1.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "2     1.0       5.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "3     1.0       5.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "4     1.0       6.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "\n",
       "   mergeMEPSNHIS  smoken   totinc  inscov  _merge  \n",
       "0            1.0     0.0  30000.0     1.0       3  \n",
       "1            1.0     0.0  25200.0     1.0       3  \n",
       "2            1.0     NaN      0.0     3.0       3  \n",
       "3            1.0     NaN      0.0     3.0       3  \n",
       "4            1.0     NaN      0.0     3.0       3  \n",
       "\n",
       "[5 rows x 125 columns]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_stata('../data_sources/meps/MEPS0004_agghealth.dta', convert_categoricals=False)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['duid',\n",
       " 'dupersid',\n",
       " 'famidyr',\n",
       " 'famszeyr',\n",
       " 'famrfpyr',\n",
       " 'region53',\n",
       " 'age',\n",
       " 'racex',\n",
       " 'hispan',\n",
       " 'marry53x',\n",
       " 'educyear',\n",
       " 'totded04',\n",
       " 'rthlth53',\n",
       " 'mnhlth53',\n",
       " 'iadlhp53',\n",
       " 'adlhlp53',\n",
       " 'bmindx53',\n",
       " 'diabdx53',\n",
       " 'asthdx53',\n",
       " 'hibpdx53',\n",
       " 'chddx53',\n",
       " 'angidx53',\n",
       " 'midx53',\n",
       " 'ohrtdx53',\n",
       " 'strkdx53',\n",
       " 'emphdx53',\n",
       " 'adsmok42',\n",
       " 'tottch04',\n",
       " 'totexp',\n",
       " 'totslf',\n",
       " 'totmcr',\n",
       " 'totmcd',\n",
       " 'totprv',\n",
       " 'totva',\n",
       " 'tottri04',\n",
       " 'totofd',\n",
       " 'totstl04',\n",
       " 'totwcp',\n",
       " 'totopr',\n",
       " 'totopu',\n",
       " 'totosr',\n",
       " 'doctim',\n",
       " 'hsptim',\n",
       " 'hspnit',\n",
       " 'perwt',\n",
       " 'famwt',\n",
       " 'male',\n",
       " 'mepsexp',\n",
       " 'mepsslf',\n",
       " 'mepsmcr',\n",
       " 'mepsmcd',\n",
       " 'mepsprv',\n",
       " 'mepsva',\n",
       " 'mepsofd',\n",
       " 'mepswcp',\n",
       " 'mepsopr',\n",
       " 'mepsopu',\n",
       " 'mepsosr',\n",
       " 'yr',\n",
       " 'totded03',\n",
       " 'tottch03',\n",
       " 'tottri03',\n",
       " 'totstl03',\n",
       " 'totded02',\n",
       " 'tottch02',\n",
       " 'tottri02',\n",
       " 'totstl02',\n",
       " 'totded01',\n",
       " 'tottch01',\n",
       " 'tottri01',\n",
       " 'totstl01',\n",
       " 'totded00',\n",
       " 'tottch00',\n",
       " 'tottri00',\n",
       " 'totstl00',\n",
       " 'cccodex',\n",
       " 'ncccodex',\n",
       " 'diabecr',\n",
       " 'hibpecr',\n",
       " 'strokecr',\n",
       " 'heartecr',\n",
       " 'lungecr',\n",
       " 'cancrecr',\n",
       " 'hearte',\n",
       " 'diabe',\n",
       " 'hibpe',\n",
       " 'stroke',\n",
       " 'lunge',\n",
       " 'adl1p',\n",
       " 'iadl1',\n",
       " 'age5659',\n",
       " 'age6064',\n",
       " 'black',\n",
       " 'educ',\n",
       " 'hsdrop',\n",
       " 'somecol',\n",
       " 'collgrad',\n",
       " 'regnth',\n",
       " 'regmid',\n",
       " 'regwst',\n",
       " 'widowed',\n",
       " 'single',\n",
       " 'bmi',\n",
       " 'obese',\n",
       " 'overwt',\n",
       " 'normalwt',\n",
       " 'underwt',\n",
       " 'wtstate',\n",
       " 'nh_hhx',\n",
       " 'nh_px',\n",
       " 'panel',\n",
       " 'nhisyr',\n",
       " 'nh_hibpe',\n",
       " 'nh_hearte',\n",
       " 'nh_stroke',\n",
       " 'nh_diabe',\n",
       " 'nh_lunge',\n",
       " 'nh_cancre',\n",
       " 'smokev',\n",
       " 'nh_smoken',\n",
       " 'mergeMEPSNHIS',\n",
       " 'smoken',\n",
       " 'totinc',\n",
       " 'inscov',\n",
       " '_merge']"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.columns.to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[df['yr'].isin([2003,2004,2005,2006,2007])]\n",
    "df = df[(df.age>=35)&(df.age<=75)]\n",
    "df = df[df.totinc>10e3]\n",
    "df = df[~df.mepsexp.isna()]\n",
    "data = df[['totinc','adl1p','age','mepsexp','totmcd','yr','male','obese','diabe','hearte','hibpe','stroke','lunge','perwt','smoken','smokev']].dropna(axis=0)\n",
    "data['loginc'] = np.log(data['totinc'])\n",
    "data['m'] = data['mepsexp'] \n",
    "data['logm'] = np.log(data.loc[data.m>0,'m'].min() + data['m'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "42.62466796875"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.01*data['m'].mean() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/users/loulou/anaconda3/lib/python3.8/site-packages/statsmodels/base/model.py:127: ValueWarning: unknown kwargs ['var_weights']\n",
      "  warnings.warn(msg, ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "model = wls('loginc ~ C(age) + C(yr) ',data=data,var_weights=np.asarray(data['perwt'])).fit()\n",
    "# + smoken + smokev + obese + diabe + stroke + hibpe + lunge + hearte + male \n",
    "data.loc[:,'ey'] = model.resid.to_list()\n",
    "data.ey = (data.ey - data.ey.mean())/data.ey.std()\n",
    "data['ey'] = data['ey'].clip(lower=-2.5,upper=2.5)\n",
    "data['ey10'] = pd.cut(data['ey'],bins=10,labels=[x for x in range(1,11)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "for q in range(1,11):\n",
    "\tdata['ey10_'+str(q)] = np.where(data['ey10']==q,1,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>totinc</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>39530.646980</td>\n",
       "      <td>29307.996623</td>\n",
       "      <td>10001.000000</td>\n",
       "      <td>19410.500000</td>\n",
       "      <td>30834.500000</td>\n",
       "      <td>50118.750000</td>\n",
       "      <td>279668.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>adl1p</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.011678</td>\n",
       "      <td>0.107438</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>51.375365</td>\n",
       "      <td>11.167422</td>\n",
       "      <td>35.000000</td>\n",
       "      <td>42.000000</td>\n",
       "      <td>50.000000</td>\n",
       "      <td>60.000000</td>\n",
       "      <td>75.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mepsexp</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>4262.466797</td>\n",
       "      <td>9039.359375</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>361.830154</td>\n",
       "      <td>1484.150269</td>\n",
       "      <td>4296.375122</td>\n",
       "      <td>205658.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>totmcd</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>260.089939</td>\n",
       "      <td>2611.798043</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>113045.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>yr</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>2003.534302</td>\n",
       "      <td>0.500991</td>\n",
       "      <td>2003.000000</td>\n",
       "      <td>2003.000000</td>\n",
       "      <td>2004.000000</td>\n",
       "      <td>2004.000000</td>\n",
       "      <td>2004.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>male</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.458567</td>\n",
       "      <td>0.498305</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>obese</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.299349</td>\n",
       "      <td>0.457999</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diabe</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.084101</td>\n",
       "      <td>0.277563</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hearte</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.104424</td>\n",
       "      <td>0.305836</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hibpe</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.310128</td>\n",
       "      <td>0.462586</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>stroke</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.024590</td>\n",
       "      <td>0.154888</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>lunge</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.017965</td>\n",
       "      <td>0.132829</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>perwt</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>10825.980579</td>\n",
       "      <td>5811.769553</td>\n",
       "      <td>711.079529</td>\n",
       "      <td>6759.226783</td>\n",
       "      <td>10100.714239</td>\n",
       "      <td>13874.193993</td>\n",
       "      <td>57180.156250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>smoken</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.222883</td>\n",
       "      <td>0.416204</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>smokev</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.482596</td>\n",
       "      <td>0.499732</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>loginc</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>10.368490</td>\n",
       "      <td>0.642292</td>\n",
       "      <td>9.210440</td>\n",
       "      <td>9.873569</td>\n",
       "      <td>10.336389</td>\n",
       "      <td>10.822150</td>\n",
       "      <td>12.541358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>m</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>4262.466797</td>\n",
       "      <td>9039.359375</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>361.830154</td>\n",
       "      <td>1484.150269</td>\n",
       "      <td>4296.375122</td>\n",
       "      <td>205658.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>logm</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>6.714816</td>\n",
       "      <td>2.538181</td>\n",
       "      <td>0.734520</td>\n",
       "      <td>5.896919</td>\n",
       "      <td>7.304001</td>\n",
       "      <td>8.366012</td>\n",
       "      <td>12.233980</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>-0.001244</td>\n",
       "      <td>0.996690</td>\n",
       "      <td>-2.028125</td>\n",
       "      <td>-0.764337</td>\n",
       "      <td>-0.054477</td>\n",
       "      <td>0.699739</td>\n",
       "      <td>2.500000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_1</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.045587</td>\n",
       "      <td>0.208600</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_2</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.098585</td>\n",
       "      <td>0.298121</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_3</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.137660</td>\n",
       "      <td>0.344562</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_4</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.159780</td>\n",
       "      <td>0.366422</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_5</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.163822</td>\n",
       "      <td>0.370135</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_6</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.141815</td>\n",
       "      <td>0.348879</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_7</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.117224</td>\n",
       "      <td>0.321705</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_8</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.075118</td>\n",
       "      <td>0.263596</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_9</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.028520</td>\n",
       "      <td>0.166463</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_10</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.031889</td>\n",
       "      <td>0.175713</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          count          mean           std           min           25%  \\\n",
       "totinc   8906.0  39530.646980  29307.996623  10001.000000  19410.500000   \n",
       "adl1p    8906.0      0.011678      0.107438      0.000000      0.000000   \n",
       "age      8906.0     51.375365     11.167422     35.000000     42.000000   \n",
       "mepsexp  8906.0   4262.466797   9039.359375      0.000000    361.830154   \n",
       "totmcd   8906.0    260.089939   2611.798043      0.000000      0.000000   \n",
       "yr       8906.0   2003.534302      0.500991   2003.000000   2003.000000   \n",
       "male     8906.0      0.458567      0.498305      0.000000      0.000000   \n",
       "obese    8906.0      0.299349      0.457999      0.000000      0.000000   \n",
       "diabe    8906.0      0.084101      0.277563      0.000000      0.000000   \n",
       "hearte   8906.0      0.104424      0.305836      0.000000      0.000000   \n",
       "hibpe    8906.0      0.310128      0.462586      0.000000      0.000000   \n",
       "stroke   8906.0      0.024590      0.154888      0.000000      0.000000   \n",
       "lunge    8906.0      0.017965      0.132829      0.000000      0.000000   \n",
       "perwt    8906.0  10825.980579   5811.769553    711.079529   6759.226783   \n",
       "smoken   8906.0      0.222883      0.416204      0.000000      0.000000   \n",
       "smokev   8906.0      0.482596      0.499732      0.000000      0.000000   \n",
       "loginc   8906.0     10.368490      0.642292      9.210440      9.873569   \n",
       "m        8906.0   4262.466797   9039.359375      0.000000    361.830154   \n",
       "logm     8906.0      6.714816      2.538181      0.734520      5.896919   \n",
       "ey       8906.0     -0.001244      0.996690     -2.028125     -0.764337   \n",
       "ey10_1   8906.0      0.045587      0.208600      0.000000      0.000000   \n",
       "ey10_2   8906.0      0.098585      0.298121      0.000000      0.000000   \n",
       "ey10_3   8906.0      0.137660      0.344562      0.000000      0.000000   \n",
       "ey10_4   8906.0      0.159780      0.366422      0.000000      0.000000   \n",
       "ey10_5   8906.0      0.163822      0.370135      0.000000      0.000000   \n",
       "ey10_6   8906.0      0.141815      0.348879      0.000000      0.000000   \n",
       "ey10_7   8906.0      0.117224      0.321705      0.000000      0.000000   \n",
       "ey10_8   8906.0      0.075118      0.263596      0.000000      0.000000   \n",
       "ey10_9   8906.0      0.028520      0.166463      0.000000      0.000000   \n",
       "ey10_10  8906.0      0.031889      0.175713      0.000000      0.000000   \n",
       "\n",
       "                  50%           75%            max  \n",
       "totinc   30834.500000  50118.750000  279668.000000  \n",
       "adl1p        0.000000      0.000000       1.000000  \n",
       "age         50.000000     60.000000      75.000000  \n",
       "mepsexp   1484.150269   4296.375122  205658.000000  \n",
       "totmcd       0.000000      0.000000  113045.000000  \n",
       "yr        2004.000000   2004.000000    2004.000000  \n",
       "male         0.000000      1.000000       1.000000  \n",
       "obese        0.000000      1.000000       1.000000  \n",
       "diabe        0.000000      0.000000       1.000000  \n",
       "hearte       0.000000      0.000000       1.000000  \n",
       "hibpe        0.000000      1.000000       1.000000  \n",
       "stroke       0.000000      0.000000       1.000000  \n",
       "lunge        0.000000      0.000000       1.000000  \n",
       "perwt    10100.714239  13874.193993   57180.156250  \n",
       "smoken       0.000000      0.000000       1.000000  \n",
       "smokev       0.000000      1.000000       1.000000  \n",
       "loginc      10.336389     10.822150      12.541358  \n",
       "m         1484.150269   4296.375122  205658.000000  \n",
       "logm         7.304001      8.366012      12.233980  \n",
       "ey          -0.054477      0.699739       2.500000  \n",
       "ey10_1       0.000000      0.000000       1.000000  \n",
       "ey10_2       0.000000      0.000000       1.000000  \n",
       "ey10_3       0.000000      0.000000       1.000000  \n",
       "ey10_4       0.000000      0.000000       1.000000  \n",
       "ey10_5       0.000000      0.000000       1.000000  \n",
       "ey10_6       0.000000      0.000000       1.000000  \n",
       "ey10_7       0.000000      0.000000       1.000000  \n",
       "ey10_8       0.000000      0.000000       1.000000  \n",
       "ey10_9       0.000000      0.000000       1.000000  \n",
       "ey10_10      0.000000      0.000000       1.000000  "
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.describe().transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/users/loulou/anaconda3/lib/python3.8/site-packages/statsmodels/base/model.py:127: ValueWarning: unknown kwargs ['var_weights']\n",
      "  warnings.warn(msg, ValueWarning)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>WLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>logm</td>       <th>  R-squared:         </th> <td>   0.113</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>WLS</td>       <th>  Adj. R-squared:    </th> <td>   0.108</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   22.53</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 27 May 2022</td> <th>  Prob (F-statistic):</th> <td>5.21e-190</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>09:09:47</td>     <th>  Log-Likelihood:    </th> <td> -20399.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  8906</td>      <th>  AIC:               </th> <td>4.090e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  8855</td>      <th>  BIC:               </th> <td>4.126e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    50</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    5.7718</td> <td>    0.158</td> <td>   36.427</td> <td> 0.000</td> <td>    5.461</td> <td>    6.082</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.36.0]</th>  <td>    0.0733</td> <td>    0.200</td> <td>    0.366</td> <td> 0.714</td> <td>   -0.319</td> <td>    0.466</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.37.0]</th>  <td>   -0.1080</td> <td>    0.202</td> <td>   -0.535</td> <td> 0.593</td> <td>   -0.504</td> <td>    0.288</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.38.0]</th>  <td>    0.0320</td> <td>    0.201</td> <td>    0.159</td> <td> 0.873</td> <td>   -0.362</td> <td>    0.426</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.39.0]</th>  <td>   -0.0883</td> <td>    0.201</td> <td>   -0.439</td> <td> 0.661</td> <td>   -0.483</td> <td>    0.306</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.40.0]</th>  <td>    0.2469</td> <td>    0.199</td> <td>    1.239</td> <td> 0.215</td> <td>   -0.144</td> <td>    0.637</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.41.0]</th>  <td>    0.2374</td> <td>    0.195</td> <td>    1.215</td> <td> 0.224</td> <td>   -0.146</td> <td>    0.620</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.42.0]</th>  <td>    0.2928</td> <td>    0.203</td> <td>    1.444</td> <td> 0.149</td> <td>   -0.105</td> <td>    0.690</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.43.0]</th>  <td>    0.2471</td> <td>    0.204</td> <td>    1.209</td> <td> 0.227</td> <td>   -0.153</td> <td>    0.648</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.44.0]</th>  <td>    0.5108</td> <td>    0.203</td> <td>    2.520</td> <td> 0.012</td> <td>    0.114</td> <td>    0.908</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.45.0]</th>  <td>    0.8884</td> <td>    0.200</td> <td>    4.450</td> <td> 0.000</td> <td>    0.497</td> <td>    1.280</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.46.0]</th>  <td>    0.8845</td> <td>    0.208</td> <td>    4.244</td> <td> 0.000</td> <td>    0.476</td> <td>    1.293</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.47.0]</th>  <td>    0.6933</td> <td>    0.203</td> <td>    3.414</td> <td> 0.001</td> <td>    0.295</td> <td>    1.091</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.48.0]</th>  <td>    0.7154</td> <td>    0.208</td> <td>    3.444</td> <td> 0.001</td> <td>    0.308</td> <td>    1.123</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.49.0]</th>  <td>    0.8715</td> <td>    0.208</td> <td>    4.191</td> <td> 0.000</td> <td>    0.464</td> <td>    1.279</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.50.0]</th>  <td>    0.9586</td> <td>    0.207</td> <td>    4.638</td> <td> 0.000</td> <td>    0.553</td> <td>    1.364</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.51.0]</th>  <td>    1.1212</td> <td>    0.211</td> <td>    5.302</td> <td> 0.000</td> <td>    0.707</td> <td>    1.536</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.52.0]</th>  <td>    1.1587</td> <td>    0.208</td> <td>    5.576</td> <td> 0.000</td> <td>    0.751</td> <td>    1.566</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.53.0]</th>  <td>    1.1226</td> <td>    0.215</td> <td>    5.211</td> <td> 0.000</td> <td>    0.700</td> <td>    1.545</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.54.0]</th>  <td>    1.3132</td> <td>    0.213</td> <td>    6.154</td> <td> 0.000</td> <td>    0.895</td> <td>    1.731</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.55.0]</th>  <td>    1.4092</td> <td>    0.217</td> <td>    6.491</td> <td> 0.000</td> <td>    0.984</td> <td>    1.835</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.56.0]</th>  <td>    1.1776</td> <td>    0.214</td> <td>    5.509</td> <td> 0.000</td> <td>    0.759</td> <td>    1.597</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.57.0]</th>  <td>    1.5739</td> <td>    0.218</td> <td>    7.232</td> <td> 0.000</td> <td>    1.147</td> <td>    2.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.58.0]</th>  <td>    1.4938</td> <td>    0.227</td> <td>    6.578</td> <td> 0.000</td> <td>    1.049</td> <td>    1.939</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.59.0]</th>  <td>    1.9002</td> <td>    0.239</td> <td>    7.956</td> <td> 0.000</td> <td>    1.432</td> <td>    2.368</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.60.0]</th>  <td>    1.6420</td> <td>    0.235</td> <td>    6.994</td> <td> 0.000</td> <td>    1.182</td> <td>    2.102</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.61.0]</th>  <td>    1.8924</td> <td>    0.236</td> <td>    8.011</td> <td> 0.000</td> <td>    1.429</td> <td>    2.355</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.62.0]</th>  <td>    2.0099</td> <td>    0.233</td> <td>    8.636</td> <td> 0.000</td> <td>    1.554</td> <td>    2.466</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.63.0]</th>  <td>    1.8638</td> <td>    0.243</td> <td>    7.661</td> <td> 0.000</td> <td>    1.387</td> <td>    2.341</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.64.0]</th>  <td>    1.8437</td> <td>    0.247</td> <td>    7.467</td> <td> 0.000</td> <td>    1.360</td> <td>    2.328</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.65.0]</th>  <td>    2.3343</td> <td>    0.248</td> <td>    9.413</td> <td> 0.000</td> <td>    1.848</td> <td>    2.820</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.66.0]</th>  <td>    1.9793</td> <td>    0.246</td> <td>    8.050</td> <td> 0.000</td> <td>    1.497</td> <td>    2.461</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.67.0]</th>  <td>    2.0425</td> <td>    0.257</td> <td>    7.935</td> <td> 0.000</td> <td>    1.538</td> <td>    2.547</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.68.0]</th>  <td>    2.2760</td> <td>    0.248</td> <td>    9.195</td> <td> 0.000</td> <td>    1.791</td> <td>    2.761</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.69.0]</th>  <td>    2.2431</td> <td>    0.259</td> <td>    8.644</td> <td> 0.000</td> <td>    1.734</td> <td>    2.752</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.70.0]</th>  <td>    2.5740</td> <td>    0.262</td> <td>    9.836</td> <td> 0.000</td> <td>    2.061</td> <td>    3.087</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.71.0]</th>  <td>    2.0607</td> <td>    0.268</td> <td>    7.692</td> <td> 0.000</td> <td>    1.536</td> <td>    2.586</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.72.0]</th>  <td>    2.6001</td> <td>    0.249</td> <td>   10.429</td> <td> 0.000</td> <td>    2.111</td> <td>    3.089</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.73.0]</th>  <td>    2.3910</td> <td>    0.250</td> <td>    9.567</td> <td> 0.000</td> <td>    1.901</td> <td>    2.881</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.74.0]</th>  <td>    2.6128</td> <td>    0.267</td> <td>    9.777</td> <td> 0.000</td> <td>    2.089</td> <td>    3.137</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.75.0]</th>  <td>    2.6292</td> <td>    0.258</td> <td>   10.182</td> <td> 0.000</td> <td>    2.123</td> <td>    3.135</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(yr)[T.2004.0]</th> <td>   -0.0472</td> <td>    0.051</td> <td>   -0.927</td> <td> 0.354</td> <td>   -0.147</td> <td>    0.053</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_1</th>          <td>   -0.3109</td> <td>    0.135</td> <td>   -2.295</td> <td> 0.022</td> <td>   -0.576</td> <td>   -0.045</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_2</th>          <td>   -0.5245</td> <td>    0.103</td> <td>   -5.110</td> <td> 0.000</td> <td>   -0.726</td> <td>   -0.323</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_3</th>          <td>   -0.4295</td> <td>    0.093</td> <td>   -4.609</td> <td> 0.000</td> <td>   -0.612</td> <td>   -0.247</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_4</th>          <td>   -0.2186</td> <td>    0.090</td> <td>   -2.440</td> <td> 0.015</td> <td>   -0.394</td> <td>   -0.043</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_6</th>          <td>    0.1822</td> <td>    0.092</td> <td>    1.973</td> <td> 0.048</td> <td>    0.001</td> <td>    0.363</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_7</th>          <td>    0.1934</td> <td>    0.097</td> <td>    1.986</td> <td> 0.047</td> <td>    0.002</td> <td>    0.384</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_8</th>          <td>    0.2131</td> <td>    0.112</td> <td>    1.900</td> <td> 0.057</td> <td>   -0.007</td> <td>    0.433</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_9</th>          <td>    0.1363</td> <td>    0.163</td> <td>    0.834</td> <td> 0.404</td> <td>   -0.184</td> <td>    0.456</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_10</th>         <td>    0.1874</td> <td>    0.156</td> <td>    1.202</td> <td> 0.229</td> <td>   -0.118</td> <td>    0.493</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>1307.653</td> <th>  Durbin-Watson:     </th> <td>   1.767</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1979.866</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td>-1.065</td>  <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 3.892</td>  <th>  Cond. No.          </th> <td>    44.8</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            WLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   logm   R-squared:                       0.113\n",
       "Model:                            WLS   Adj. R-squared:                  0.108\n",
       "Method:                 Least Squares   F-statistic:                     22.53\n",
       "Date:                Fri, 27 May 2022   Prob (F-statistic):          5.21e-190\n",
       "Time:                        09:09:47   Log-Likelihood:                -20399.\n",
       "No. Observations:                8906   AIC:                         4.090e+04\n",
       "Df Residuals:                    8855   BIC:                         4.126e+04\n",
       "Df Model:                          50                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===================================================================================\n",
       "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
       "-----------------------------------------------------------------------------------\n",
       "Intercept           5.7718      0.158     36.427      0.000       5.461       6.082\n",
       "C(age)[T.36.0]      0.0733      0.200      0.366      0.714      -0.319       0.466\n",
       "C(age)[T.37.0]     -0.1080      0.202     -0.535      0.593      -0.504       0.288\n",
       "C(age)[T.38.0]      0.0320      0.201      0.159      0.873      -0.362       0.426\n",
       "C(age)[T.39.0]     -0.0883      0.201     -0.439      0.661      -0.483       0.306\n",
       "C(age)[T.40.0]      0.2469      0.199      1.239      0.215      -0.144       0.637\n",
       "C(age)[T.41.0]      0.2374      0.195      1.215      0.224      -0.146       0.620\n",
       "C(age)[T.42.0]      0.2928      0.203      1.444      0.149      -0.105       0.690\n",
       "C(age)[T.43.0]      0.2471      0.204      1.209      0.227      -0.153       0.648\n",
       "C(age)[T.44.0]      0.5108      0.203      2.520      0.012       0.114       0.908\n",
       "C(age)[T.45.0]      0.8884      0.200      4.450      0.000       0.497       1.280\n",
       "C(age)[T.46.0]      0.8845      0.208      4.244      0.000       0.476       1.293\n",
       "C(age)[T.47.0]      0.6933      0.203      3.414      0.001       0.295       1.091\n",
       "C(age)[T.48.0]      0.7154      0.208      3.444      0.001       0.308       1.123\n",
       "C(age)[T.49.0]      0.8715      0.208      4.191      0.000       0.464       1.279\n",
       "C(age)[T.50.0]      0.9586      0.207      4.638      0.000       0.553       1.364\n",
       "C(age)[T.51.0]      1.1212      0.211      5.302      0.000       0.707       1.536\n",
       "C(age)[T.52.0]      1.1587      0.208      5.576      0.000       0.751       1.566\n",
       "C(age)[T.53.0]      1.1226      0.215      5.211      0.000       0.700       1.545\n",
       "C(age)[T.54.0]      1.3132      0.213      6.154      0.000       0.895       1.731\n",
       "C(age)[T.55.0]      1.4092      0.217      6.491      0.000       0.984       1.835\n",
       "C(age)[T.56.0]      1.1776      0.214      5.509      0.000       0.759       1.597\n",
       "C(age)[T.57.0]      1.5739      0.218      7.232      0.000       1.147       2.000\n",
       "C(age)[T.58.0]      1.4938      0.227      6.578      0.000       1.049       1.939\n",
       "C(age)[T.59.0]      1.9002      0.239      7.956      0.000       1.432       2.368\n",
       "C(age)[T.60.0]      1.6420      0.235      6.994      0.000       1.182       2.102\n",
       "C(age)[T.61.0]      1.8924      0.236      8.011      0.000       1.429       2.355\n",
       "C(age)[T.62.0]      2.0099      0.233      8.636      0.000       1.554       2.466\n",
       "C(age)[T.63.0]      1.8638      0.243      7.661      0.000       1.387       2.341\n",
       "C(age)[T.64.0]      1.8437      0.247      7.467      0.000       1.360       2.328\n",
       "C(age)[T.65.0]      2.3343      0.248      9.413      0.000       1.848       2.820\n",
       "C(age)[T.66.0]      1.9793      0.246      8.050      0.000       1.497       2.461\n",
       "C(age)[T.67.0]      2.0425      0.257      7.935      0.000       1.538       2.547\n",
       "C(age)[T.68.0]      2.2760      0.248      9.195      0.000       1.791       2.761\n",
       "C(age)[T.69.0]      2.2431      0.259      8.644      0.000       1.734       2.752\n",
       "C(age)[T.70.0]      2.5740      0.262      9.836      0.000       2.061       3.087\n",
       "C(age)[T.71.0]      2.0607      0.268      7.692      0.000       1.536       2.586\n",
       "C(age)[T.72.0]      2.6001      0.249     10.429      0.000       2.111       3.089\n",
       "C(age)[T.73.0]      2.3910      0.250      9.567      0.000       1.901       2.881\n",
       "C(age)[T.74.0]      2.6128      0.267      9.777      0.000       2.089       3.137\n",
       "C(age)[T.75.0]      2.6292      0.258     10.182      0.000       2.123       3.135\n",
       "C(yr)[T.2004.0]    -0.0472      0.051     -0.927      0.354      -0.147       0.053\n",
       "ey10_1             -0.3109      0.135     -2.295      0.022      -0.576      -0.045\n",
       "ey10_2             -0.5245      0.103     -5.110      0.000      -0.726      -0.323\n",
       "ey10_3             -0.4295      0.093     -4.609      0.000      -0.612      -0.247\n",
       "ey10_4             -0.2186      0.090     -2.440      0.015      -0.394      -0.043\n",
       "ey10_6              0.1822      0.092      1.973      0.048       0.001       0.363\n",
       "ey10_7              0.1934      0.097      1.986      0.047       0.002       0.384\n",
       "ey10_8              0.2131      0.112      1.900      0.057      -0.007       0.433\n",
       "ey10_9              0.1363      0.163      0.834      0.404      -0.184       0.456\n",
       "ey10_10             0.1874      0.156      1.202      0.229      -0.118       0.493\n",
       "==============================================================================\n",
       "Omnibus:                     1307.653   Durbin-Watson:                   1.767\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1979.866\n",
       "Skew:                          -1.065   Prob(JB):                         0.00\n",
       "Kurtosis:                       3.892   Cond. No.                         44.8\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = wls('logm ~ ey10_1 +ey10_2 + ey10_3 + ey10_4 + ey10_6 + ey10_7 + ey10_8 + ey10_9 + ey10_10 + C(age) + C(yr) ',data=data,var_weights=np.asarray(data['perwt'])).fit()\n",
    "#+ smoken + smokev + obese + diabe + stroke + hibpe + lunge + hearte + male\n",
    "model.summary()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1     0.156416\n",
       "2     0.114373\n",
       "3     0.095196\n",
       "4     0.090118\n",
       "5     0.000000\n",
       "6     0.085721\n",
       "7     0.090745\n",
       "8     0.100431\n",
       "9     0.127536\n",
       "10    0.131535\n",
       "dtype: float64"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_data_se = model.HC0_se[['ey10_'+str(q) for q in range(1,11) if q!=5]]\n",
    "grad_m_data_se.index = [x for x in range(1,11) if x!=5]\n",
    "grad_m_data_se[5] = 0\n",
    "grad_m_data_se = grad_m_data_se.sort_index()\n",
    "grad_m_data_se"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    -0.310938\n",
       "2    -0.524490\n",
       "3    -0.429465\n",
       "4    -0.218555\n",
       "5     0.000000\n",
       "6     0.182247\n",
       "7     0.193383\n",
       "8     0.213062\n",
       "9     0.136299\n",
       "10    0.187443\n",
       "dtype: float64"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_data = model.params[['ey10_'+str(q) for q in range(1,11) if q!=5]]\n",
    "grad_m_data.index = [x for x in range(1,11) if x!=5]\n",
    "grad_m_data[5] = 0\n",
    "grad_m_data = grad_m_data.sort_index()\n",
    "grad_m_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/users/loulou/anaconda3/lib/python3.8/site-packages/statsmodels/base/model.py:127: ValueWarning: unknown kwargs ['var_weights']\n",
      "  warnings.warn(msg, ValueWarning)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>WLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>logm</td>       <th>  R-squared:         </th> <td>   0.226</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>WLS</td>       <th>  Adj. R-squared:    </th> <td>   0.220</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   43.67</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 27 May 2022</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>09:09:48</td>     <th>  Log-Likelihood:    </th> <td> -19794.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  8906</td>      <th>  AIC:               </th> <td>3.971e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  8846</td>      <th>  BIC:               </th> <td>4.013e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    59</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    6.0258</td> <td>    0.152</td> <td>   39.556</td> <td> 0.000</td> <td>    5.727</td> <td>    6.324</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.36.0]</th>  <td>    0.0392</td> <td>    0.187</td> <td>    0.209</td> <td> 0.834</td> <td>   -0.328</td> <td>    0.406</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.37.0]</th>  <td>   -0.0815</td> <td>    0.189</td> <td>   -0.432</td> <td> 0.666</td> <td>   -0.452</td> <td>    0.289</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.38.0]</th>  <td>    0.0023</td> <td>    0.188</td> <td>    0.012</td> <td> 0.990</td> <td>   -0.366</td> <td>    0.371</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.39.0]</th>  <td>   -0.1169</td> <td>    0.188</td> <td>   -0.621</td> <td> 0.534</td> <td>   -0.486</td> <td>    0.252</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.40.0]</th>  <td>    0.2486</td> <td>    0.186</td> <td>    1.335</td> <td> 0.182</td> <td>   -0.116</td> <td>    0.614</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.41.0]</th>  <td>    0.1654</td> <td>    0.183</td> <td>    0.905</td> <td> 0.365</td> <td>   -0.193</td> <td>    0.524</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.42.0]</th>  <td>    0.2213</td> <td>    0.190</td> <td>    1.167</td> <td> 0.243</td> <td>   -0.151</td> <td>    0.593</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.43.0]</th>  <td>    0.1910</td> <td>    0.191</td> <td>    0.999</td> <td> 0.318</td> <td>   -0.184</td> <td>    0.566</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.44.0]</th>  <td>    0.3976</td> <td>    0.190</td> <td>    2.097</td> <td> 0.036</td> <td>    0.026</td> <td>    0.769</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.45.0]</th>  <td>    0.7522</td> <td>    0.187</td> <td>    4.027</td> <td> 0.000</td> <td>    0.386</td> <td>    1.118</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.46.0]</th>  <td>    0.6286</td> <td>    0.195</td> <td>    3.224</td> <td> 0.001</td> <td>    0.246</td> <td>    1.011</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.47.0]</th>  <td>    0.4826</td> <td>    0.190</td> <td>    2.540</td> <td> 0.011</td> <td>    0.110</td> <td>    0.855</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.48.0]</th>  <td>    0.5367</td> <td>    0.194</td> <td>    2.761</td> <td> 0.006</td> <td>    0.156</td> <td>    0.918</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.49.0]</th>  <td>    0.5856</td> <td>    0.195</td> <td>    3.008</td> <td> 0.003</td> <td>    0.204</td> <td>    0.967</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.50.0]</th>  <td>    0.5697</td> <td>    0.194</td> <td>    2.943</td> <td> 0.003</td> <td>    0.190</td> <td>    0.949</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.51.0]</th>  <td>    0.7385</td> <td>    0.198</td> <td>    3.727</td> <td> 0.000</td> <td>    0.350</td> <td>    1.127</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.52.0]</th>  <td>    0.6739</td> <td>    0.195</td> <td>    3.459</td> <td> 0.001</td> <td>    0.292</td> <td>    1.056</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.53.0]</th>  <td>    0.6313</td> <td>    0.202</td> <td>    3.125</td> <td> 0.002</td> <td>    0.235</td> <td>    1.027</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.54.0]</th>  <td>    0.7670</td> <td>    0.200</td> <td>    3.830</td> <td> 0.000</td> <td>    0.374</td> <td>    1.160</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.55.0]</th>  <td>    0.8264</td> <td>    0.204</td> <td>    4.055</td> <td> 0.000</td> <td>    0.427</td> <td>    1.226</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.56.0]</th>  <td>    0.7326</td> <td>    0.201</td> <td>    3.653</td> <td> 0.000</td> <td>    0.339</td> <td>    1.126</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.57.0]</th>  <td>    1.0893</td> <td>    0.204</td> <td>    5.329</td> <td> 0.000</td> <td>    0.689</td> <td>    1.490</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.58.0]</th>  <td>    0.9662</td> <td>    0.213</td> <td>    4.533</td> <td> 0.000</td> <td>    0.548</td> <td>    1.384</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.59.0]</th>  <td>    1.2559</td> <td>    0.224</td> <td>    5.596</td> <td> 0.000</td> <td>    0.816</td> <td>    1.696</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.60.0]</th>  <td>    0.8639</td> <td>    0.221</td> <td>    3.908</td> <td> 0.000</td> <td>    0.431</td> <td>    1.297</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.61.0]</th>  <td>    1.0555</td> <td>    0.223</td> <td>    4.741</td> <td> 0.000</td> <td>    0.619</td> <td>    1.492</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.62.0]</th>  <td>    1.3250</td> <td>    0.219</td> <td>    6.051</td> <td> 0.000</td> <td>    0.896</td> <td>    1.754</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.63.0]</th>  <td>    1.2049</td> <td>    0.229</td> <td>    5.263</td> <td> 0.000</td> <td>    0.756</td> <td>    1.654</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.64.0]</th>  <td>    1.2446</td> <td>    0.232</td> <td>    5.363</td> <td> 0.000</td> <td>    0.790</td> <td>    1.700</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.65.0]</th>  <td>    1.4242</td> <td>    0.234</td> <td>    6.084</td> <td> 0.000</td> <td>    0.965</td> <td>    1.883</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.66.0]</th>  <td>    0.9755</td> <td>    0.232</td> <td>    4.200</td> <td> 0.000</td> <td>    0.520</td> <td>    1.431</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.67.0]</th>  <td>    1.0702</td> <td>    0.243</td> <td>    4.399</td> <td> 0.000</td> <td>    0.593</td> <td>    1.547</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.68.0]</th>  <td>    1.2799</td> <td>    0.234</td> <td>    5.466</td> <td> 0.000</td> <td>    0.821</td> <td>    1.739</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.69.0]</th>  <td>    1.3089</td> <td>    0.245</td> <td>    5.345</td> <td> 0.000</td> <td>    0.829</td> <td>    1.789</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.70.0]</th>  <td>    1.5193</td> <td>    0.247</td> <td>    6.144</td> <td> 0.000</td> <td>    1.035</td> <td>    2.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.71.0]</th>  <td>    1.0835</td> <td>    0.253</td> <td>    4.281</td> <td> 0.000</td> <td>    0.587</td> <td>    1.580</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.72.0]</th>  <td>    1.3171</td> <td>    0.237</td> <td>    5.558</td> <td> 0.000</td> <td>    0.853</td> <td>    1.782</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.73.0]</th>  <td>    1.1630</td> <td>    0.237</td> <td>    4.901</td> <td> 0.000</td> <td>    0.698</td> <td>    1.628</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.74.0]</th>  <td>    1.4371</td> <td>    0.253</td> <td>    5.671</td> <td> 0.000</td> <td>    0.940</td> <td>    1.934</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.75.0]</th>  <td>    1.3932</td> <td>    0.245</td> <td>    5.676</td> <td> 0.000</td> <td>    0.912</td> <td>    1.874</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(yr)[T.2004.0]</th> <td>   -0.0580</td> <td>    0.048</td> <td>   -1.218</td> <td> 0.223</td> <td>   -0.151</td> <td>    0.035</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_1</th>          <td>   -0.5309</td> <td>    0.127</td> <td>   -4.177</td> <td> 0.000</td> <td>   -0.780</td> <td>   -0.282</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_2</th>          <td>   -0.6819</td> <td>    0.096</td> <td>   -7.088</td> <td> 0.000</td> <td>   -0.870</td> <td>   -0.493</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_3</th>          <td>   -0.5519</td> <td>    0.087</td> <td>   -6.327</td> <td> 0.000</td> <td>   -0.723</td> <td>   -0.381</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_4</th>          <td>   -0.2330</td> <td>    0.084</td> <td>   -2.780</td> <td> 0.005</td> <td>   -0.397</td> <td>   -0.069</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_6</th>          <td>    0.2285</td> <td>    0.086</td> <td>    2.645</td> <td> 0.008</td> <td>    0.059</td> <td>    0.398</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_7</th>          <td>    0.2879</td> <td>    0.091</td> <td>    3.156</td> <td> 0.002</td> <td>    0.109</td> <td>    0.467</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_8</th>          <td>    0.4109</td> <td>    0.105</td> <td>    3.904</td> <td> 0.000</td> <td>    0.205</td> <td>    0.617</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_9</th>          <td>    0.4310</td> <td>    0.153</td> <td>    2.813</td> <td> 0.005</td> <td>    0.131</td> <td>    0.731</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_10</th>         <td>    0.5187</td> <td>    0.147</td> <td>    3.540</td> <td> 0.000</td> <td>    0.231</td> <td>    0.806</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smoken</th>          <td>   -0.2906</td> <td>    0.068</td> <td>   -4.297</td> <td> 0.000</td> <td>   -0.423</td> <td>   -0.158</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokev</th>          <td>    0.2643</td> <td>    0.056</td> <td>    4.710</td> <td> 0.000</td> <td>    0.154</td> <td>    0.374</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>obese</th>           <td>    0.0901</td> <td>    0.054</td> <td>    1.662</td> <td> 0.097</td> <td>   -0.016</td> <td>    0.196</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>diabe</th>           <td>    0.7981</td> <td>    0.091</td> <td>    8.725</td> <td> 0.000</td> <td>    0.619</td> <td>    0.977</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>stroke</th>          <td>    0.8748</td> <td>    0.158</td> <td>    5.537</td> <td> 0.000</td> <td>    0.565</td> <td>    1.185</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hibpe</th>           <td>    1.0223</td> <td>    0.058</td> <td>   17.647</td> <td> 0.000</td> <td>    0.909</td> <td>    1.136</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>lunge</th>           <td>    0.7732</td> <td>    0.184</td> <td>    4.197</td> <td> 0.000</td> <td>    0.412</td> <td>    1.134</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hearte</th>          <td>    0.9840</td> <td>    0.083</td> <td>   11.867</td> <td> 0.000</td> <td>    0.821</td> <td>    1.147</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>male</th>            <td>   -0.9986</td> <td>    0.049</td> <td>  -20.546</td> <td> 0.000</td> <td>   -1.094</td> <td>   -0.903</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>1051.272</td> <th>  Durbin-Watson:     </th> <td>   1.805</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1484.205</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td>-0.918</td>  <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 3.792</td>  <th>  Cond. No.          </th> <td>    56.2</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            WLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   logm   R-squared:                       0.226\n",
       "Model:                            WLS   Adj. R-squared:                  0.220\n",
       "Method:                 Least Squares   F-statistic:                     43.67\n",
       "Date:                Fri, 27 May 2022   Prob (F-statistic):               0.00\n",
       "Time:                        09:09:48   Log-Likelihood:                -19794.\n",
       "No. Observations:                8906   AIC:                         3.971e+04\n",
       "Df Residuals:                    8846   BIC:                         4.013e+04\n",
       "Df Model:                          59                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===================================================================================\n",
       "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
       "-----------------------------------------------------------------------------------\n",
       "Intercept           6.0258      0.152     39.556      0.000       5.727       6.324\n",
       "C(age)[T.36.0]      0.0392      0.187      0.209      0.834      -0.328       0.406\n",
       "C(age)[T.37.0]     -0.0815      0.189     -0.432      0.666      -0.452       0.289\n",
       "C(age)[T.38.0]      0.0023      0.188      0.012      0.990      -0.366       0.371\n",
       "C(age)[T.39.0]     -0.1169      0.188     -0.621      0.534      -0.486       0.252\n",
       "C(age)[T.40.0]      0.2486      0.186      1.335      0.182      -0.116       0.614\n",
       "C(age)[T.41.0]      0.1654      0.183      0.905      0.365      -0.193       0.524\n",
       "C(age)[T.42.0]      0.2213      0.190      1.167      0.243      -0.151       0.593\n",
       "C(age)[T.43.0]      0.1910      0.191      0.999      0.318      -0.184       0.566\n",
       "C(age)[T.44.0]      0.3976      0.190      2.097      0.036       0.026       0.769\n",
       "C(age)[T.45.0]      0.7522      0.187      4.027      0.000       0.386       1.118\n",
       "C(age)[T.46.0]      0.6286      0.195      3.224      0.001       0.246       1.011\n",
       "C(age)[T.47.0]      0.4826      0.190      2.540      0.011       0.110       0.855\n",
       "C(age)[T.48.0]      0.5367      0.194      2.761      0.006       0.156       0.918\n",
       "C(age)[T.49.0]      0.5856      0.195      3.008      0.003       0.204       0.967\n",
       "C(age)[T.50.0]      0.5697      0.194      2.943      0.003       0.190       0.949\n",
       "C(age)[T.51.0]      0.7385      0.198      3.727      0.000       0.350       1.127\n",
       "C(age)[T.52.0]      0.6739      0.195      3.459      0.001       0.292       1.056\n",
       "C(age)[T.53.0]      0.6313      0.202      3.125      0.002       0.235       1.027\n",
       "C(age)[T.54.0]      0.7670      0.200      3.830      0.000       0.374       1.160\n",
       "C(age)[T.55.0]      0.8264      0.204      4.055      0.000       0.427       1.226\n",
       "C(age)[T.56.0]      0.7326      0.201      3.653      0.000       0.339       1.126\n",
       "C(age)[T.57.0]      1.0893      0.204      5.329      0.000       0.689       1.490\n",
       "C(age)[T.58.0]      0.9662      0.213      4.533      0.000       0.548       1.384\n",
       "C(age)[T.59.0]      1.2559      0.224      5.596      0.000       0.816       1.696\n",
       "C(age)[T.60.0]      0.8639      0.221      3.908      0.000       0.431       1.297\n",
       "C(age)[T.61.0]      1.0555      0.223      4.741      0.000       0.619       1.492\n",
       "C(age)[T.62.0]      1.3250      0.219      6.051      0.000       0.896       1.754\n",
       "C(age)[T.63.0]      1.2049      0.229      5.263      0.000       0.756       1.654\n",
       "C(age)[T.64.0]      1.2446      0.232      5.363      0.000       0.790       1.700\n",
       "C(age)[T.65.0]      1.4242      0.234      6.084      0.000       0.965       1.883\n",
       "C(age)[T.66.0]      0.9755      0.232      4.200      0.000       0.520       1.431\n",
       "C(age)[T.67.0]      1.0702      0.243      4.399      0.000       0.593       1.547\n",
       "C(age)[T.68.0]      1.2799      0.234      5.466      0.000       0.821       1.739\n",
       "C(age)[T.69.0]      1.3089      0.245      5.345      0.000       0.829       1.789\n",
       "C(age)[T.70.0]      1.5193      0.247      6.144      0.000       1.035       2.004\n",
       "C(age)[T.71.0]      1.0835      0.253      4.281      0.000       0.587       1.580\n",
       "C(age)[T.72.0]      1.3171      0.237      5.558      0.000       0.853       1.782\n",
       "C(age)[T.73.0]      1.1630      0.237      4.901      0.000       0.698       1.628\n",
       "C(age)[T.74.0]      1.4371      0.253      5.671      0.000       0.940       1.934\n",
       "C(age)[T.75.0]      1.3932      0.245      5.676      0.000       0.912       1.874\n",
       "C(yr)[T.2004.0]    -0.0580      0.048     -1.218      0.223      -0.151       0.035\n",
       "ey10_1             -0.5309      0.127     -4.177      0.000      -0.780      -0.282\n",
       "ey10_2             -0.6819      0.096     -7.088      0.000      -0.870      -0.493\n",
       "ey10_3             -0.5519      0.087     -6.327      0.000      -0.723      -0.381\n",
       "ey10_4             -0.2330      0.084     -2.780      0.005      -0.397      -0.069\n",
       "ey10_6              0.2285      0.086      2.645      0.008       0.059       0.398\n",
       "ey10_7              0.2879      0.091      3.156      0.002       0.109       0.467\n",
       "ey10_8              0.4109      0.105      3.904      0.000       0.205       0.617\n",
       "ey10_9              0.4310      0.153      2.813      0.005       0.131       0.731\n",
       "ey10_10             0.5187      0.147      3.540      0.000       0.231       0.806\n",
       "smoken             -0.2906      0.068     -4.297      0.000      -0.423      -0.158\n",
       "smokev              0.2643      0.056      4.710      0.000       0.154       0.374\n",
       "obese               0.0901      0.054      1.662      0.097      -0.016       0.196\n",
       "diabe               0.7981      0.091      8.725      0.000       0.619       0.977\n",
       "stroke              0.8748      0.158      5.537      0.000       0.565       1.185\n",
       "hibpe               1.0223      0.058     17.647      0.000       0.909       1.136\n",
       "lunge               0.7732      0.184      4.197      0.000       0.412       1.134\n",
       "hearte              0.9840      0.083     11.867      0.000       0.821       1.147\n",
       "male               -0.9986      0.049    -20.546      0.000      -1.094      -0.903\n",
       "==============================================================================\n",
       "Omnibus:                     1051.272   Durbin-Watson:                   1.805\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1484.205\n",
       "Skew:                          -0.918   Prob(JB):                         0.00\n",
       "Kurtosis:                       3.792   Cond. No.                         56.2\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = wls('logm ~ ey10_1 +ey10_2 + ey10_3 + ey10_4 + ey10_6 + ey10_7 + ey10_8 + ey10_9 + ey10_10 + C(age) + C(yr)  + smoken + smokev + obese + diabe + stroke + hibpe + lunge + hearte + male',data=data,var_weights=np.asarray(data['perwt'])).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    -0.530896\n",
       "2    -0.681907\n",
       "3    -0.551857\n",
       "4    -0.233032\n",
       "5     0.000000\n",
       "6     0.228472\n",
       "7     0.287888\n",
       "8     0.410882\n",
       "9     0.431045\n",
       "10    0.518710\n",
       "dtype: float64"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_data_h = model.params[['ey10_'+str(q) for q in range(1,11) if q!=5]]\n",
    "grad_m_data_h.index = [x for x in range(1,11) if x!=5]\n",
    "grad_m_data_h[5] = 0\n",
    "grad_m_data_h = grad_m_data_h.sort_index()\n",
    "grad_m_data_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = pd.DataFrame(index=[x for x in range(1,11)],columns=['sim','mean','se'])\n",
    "table['sim'] = grad_m_sim \n",
    "table['mean_h'] = grad_m_data_h\n",
    "table['mean'] = grad_m_data\n",
    "table['se'] = grad_m_data_se "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sim</th>\n",
       "      <th>mean</th>\n",
       "      <th>se</th>\n",
       "      <th>mean_h</th>\n",
       "      <th>up</th>\n",
       "      <th>low</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.739822</td>\n",
       "      <td>-0.310938</td>\n",
       "      <td>0.156416</td>\n",
       "      <td>-0.530896</td>\n",
       "      <td>-0.004362</td>\n",
       "      <td>-0.617514</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.663559</td>\n",
       "      <td>-0.524490</td>\n",
       "      <td>0.114373</td>\n",
       "      <td>-0.681907</td>\n",
       "      <td>-0.300318</td>\n",
       "      <td>-0.748661</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.559215</td>\n",
       "      <td>-0.429465</td>\n",
       "      <td>0.095196</td>\n",
       "      <td>-0.551857</td>\n",
       "      <td>-0.242880</td>\n",
       "      <td>-0.616049</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.366673</td>\n",
       "      <td>-0.218555</td>\n",
       "      <td>0.090118</td>\n",
       "      <td>-0.233032</td>\n",
       "      <td>-0.041925</td>\n",
       "      <td>-0.395186</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.055878</td>\n",
       "      <td>0.182247</td>\n",
       "      <td>0.085721</td>\n",
       "      <td>0.228472</td>\n",
       "      <td>0.350260</td>\n",
       "      <td>0.014234</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>-0.032321</td>\n",
       "      <td>0.193383</td>\n",
       "      <td>0.090745</td>\n",
       "      <td>0.287888</td>\n",
       "      <td>0.371243</td>\n",
       "      <td>0.015524</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>-0.017302</td>\n",
       "      <td>0.213062</td>\n",
       "      <td>0.100431</td>\n",
       "      <td>0.410882</td>\n",
       "      <td>0.409907</td>\n",
       "      <td>0.016217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.036369</td>\n",
       "      <td>0.136299</td>\n",
       "      <td>0.127536</td>\n",
       "      <td>0.431045</td>\n",
       "      <td>0.386269</td>\n",
       "      <td>-0.113671</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.121999</td>\n",
       "      <td>0.187443</td>\n",
       "      <td>0.131535</td>\n",
       "      <td>0.518710</td>\n",
       "      <td>0.445253</td>\n",
       "      <td>-0.070366</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         sim      mean        se    mean_h        up       low\n",
       "1  -0.739822 -0.310938  0.156416 -0.530896 -0.004362 -0.617514\n",
       "2  -0.663559 -0.524490  0.114373 -0.681907 -0.300318 -0.748661\n",
       "3  -0.559215 -0.429465  0.095196 -0.551857 -0.242880 -0.616049\n",
       "4  -0.366673 -0.218555  0.090118 -0.233032 -0.041925 -0.395186\n",
       "5   0.000000  0.000000  0.000000  0.000000  0.000000  0.000000\n",
       "6   0.055878  0.182247  0.085721  0.228472  0.350260  0.014234\n",
       "7  -0.032321  0.193383  0.090745  0.287888  0.371243  0.015524\n",
       "8  -0.017302  0.213062  0.100431  0.410882  0.409907  0.016217\n",
       "9   0.036369  0.136299  0.127536  0.431045  0.386269 -0.113671\n",
       "10  0.121999  0.187443  0.131535  0.518710  0.445253 -0.070366"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table['up'] = table['mean'] + 1.96*table['se']\n",
    "table['low'] = table['mean'] - 1.96*table['se']\n",
    "table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/h9/sn9_m3953j93s5pv8q78dfvh0000gn/T/ipykernel_6728/3838045957.py:5: UserWarning: linestyle is redundantly defined by the 'linestyle' keyword argument and the fmt string \"-D\" (-> linestyle='-'). The keyword argument will take precedence.\n",
      "  ax.plot(table.index,table['mean_h'],'-D',color='r',linestyle=':',label='Data (controls for health)')\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABq3ElEQVR4nO2dd3iUxfbHP5NNbxAICSV0gat0pCoKCHJVEFH0KiBdsaGigu13FcvVq4JYrl4URUQEQb0IiL2goAYUFFFRQCFASCWkZ0uye35/bCFlN9kku2nM53n22bfMO3N2s3m/78ycOUeJCBqNRqPReCKgvg3QaDQaTcNGC4VGo9FoKkULhUaj0WgqRQuFRqPRaCpFC4VGo9FoKiWwvg3wB7GxsdKpU6f6NkOj0WgaDbt37z4hIq3cnWuSQtGpUyd27dpV32ZoNBpNo0EpdcTTOT30pNFoNJpK0UKh0Wg0mkrRQqHRaDSaSmmScxTuKC4uJjk5GZPJVN+maE5TQkNDSUhIICgoqL5N0WiqxWkjFMnJyURFRdGpUyeUUvVtjuY0Q0TIysoiOTmZzp0717c5mqbI1q0waxasXAmjRvm06tNm6MlkMtGyZUstEpp6QSlFy5YtdY9W4x+2boXx4+HIEfv71q0+rf60EQpAi4SmXtG/P41fcIpEUZF9v6jI52JxWgmFRqPRNCnKi4QTH4uFFgqNRqNprMyaVVEknBQV2c/7AC0UHlizBjp1goAA+/uaNbWv02Aw0K9fP3r27Enfvn1ZunQpNput0muSkpJYu3ZttdsyGo2MGDECq9VaU3N9zrPPPkuRpx91JURGRlZ6fsyYMWRnZ9fULI2m8fHJJ7Bxo33iOjzcfZnwcPt5H6CFwg1r1sDcufZ5IRH7+9y5tReLsLAw9uzZw2+//cZnn33Ghx9+yMMPP1zpNTUVitdee40rrrgCg8FQU3N9TmVCURtBmzZtGv/9739rfL1G0ygwGk9t//vf8OSTdu+mLVsqikV4uP24r7yfRKTJvc4++2wpz759+1zbt98uMmKE51dIiIhdIsq+QkI8X3P77RWarEBERESZ/b/++ktatGghNptNDh8+LMOHD5f+/ftL//795dtvvxURkSFDhkh0dLT07dtXli5d6rFceYYNGyaHDx8WEZGtW7fKiBEjZNKkSdKjRw+ZMmWK2Gw2ERH5/PPPpV+/ftKrVy+ZNWuWmEymCnUdPHhQRo8eLX369JH+/fvLn3/+KTabTRYsWCA9e/aUXr16ybp16ypt67nnnpOgoCDp1auXjBw50vV9PPDAAzJ48GDZvn27PP3009KzZ0/p2bOnPPPMMxW+t5SUFDnvvPOkb9++0rNnT9m2bZuIiJw8eVJ69uxZ9R+gAVD6d6jReM0LL4i0bCmSn2/fP3pUxGw+df7LL0XCw+03qvBw+341AXaJh3tqvd/U/fGqrVC4Ewnny5dCISLSvHlzSUtLk8LCQjEajSIicuDAAXF+hq1bt8q4ceNc5T2VK43ZbJb4+HjX/tatWyU6OlqOHTsmVqtVhg4dKtu3bxej0SgJCQmyf/9+ERGZNm1amRu0k8GDB8uGDRtERMRoNEphYaG8++67MmbMGCkpKZG0tDRp3769pKSkeGxLRKRjx46SmZnpqheQ9evXi4jIrl27pFevXlJQUCD5+fly1llnyY8//ljme1uyZIn861//EhGRkpISycvLc9V1xhlnyIkTJzx/+Q0ELRQar8jOFlmyROTIEfv+jh0id94pUtlv/MsvRTp2rJFIiFQuFKfNgrvSPPts5ec7dbIPN5WnY0f46ivf2mL/+9hXjs+bN489e/ZgMBg4cOCA2/LelDtx4gTNmzcvc2zw4MEkJCQA0K9fP5KSkoiKiqJz5850794dgBkzZvDiiy8yf/5813X5+fkcP36cyy+/HLCvLgb45ptvmDx5MgaDgfj4eEaMGMEPP/xAdHS027aGDx9ewU6DwcCkSZNc9V1++eVEREQAcMUVV7B9+3b69+/vKj9o0CBmz55NcXExEydOpF+/fq5zcXFxpKSk0LJlS/dftEbT0BGBwkKIjITcXFi40L59ww0wZIj9VRmjRkFSkl9M03MUbnjsMfdDfo895tt2Dh06hMFgIC4ujmeeeYb4+Hh+/vlndu3ahcVicXuNN+XCwsIqLOwKCQlxbRsMBkpKSlwiVRmeylR2rbu23BEaGuqaQ/HGlvPPP59t27bRrl07pk2bxhtvvOE6ZzKZCAsLq7IOjaZBIgLDhsG8efb9jh3h8GG7SDQAtFC4YepUWL7c/rdSyv6+fLn9uK/IzMzkxhtvZN68eSilyM3NpU2bNgQEBLB69WrX5G5UVBT5+fmu6zyVK01MTAxWq7XKVcB/+9vfSEpK4s8//wRg9erVjBgxokyZ6OhoEhIS2LhxIwBms5mioiLOP/981q9fj9VqJTMzk23btjF48OBK2yv/WUpz/vnns3HjRoqKiigsLOS9997jvPPOK1PmyJEjxMXFcf311zNnzhx+/PFHwC4yaWlp6GRVmkbFjz/aJ6TBfqO56iq48MJT5zt2rB+73HBaDj15w9SpvhUGsLus9uvXj+LiYgIDA5k2bRp33nknADfffDOTJk3inXfeYdSoUa4hmD59+hAYGEjfvn2ZOXOmx3LlGTt2LN988w1jxozxaE9oaCgrV67kqquuoqSkhEGDBnHjjTdWKLd69WpuuOEGHnzwQYKCgnjnnXe4/PLLSUxMpG/fviileOqpp2jdujV//PGHx/bmzp3LxRdfTJs2bdhabiHQgAEDmDlzpktsrrvuujLDTgBfffUVixcvJigoiMjISFePYvfu3QwdOpTAQP1z1jRw8vMhIsLud//JJ/D44zBnDsTGwl131bd1HlHedPkbGwMHDpTyGe5+//13zjzzzHqyqO756aefWLp0KatXr65vU/zO7bffzoQJExg9enR9m1Ilp9vvUFOKnTth9GjYtMn+npdnPx4dXb92OVBK7RaRge7O6aGnJkr//v0ZNWpUg1pw5y969erVKERC00TZutXuAVM+XIbVCitW2NczAPTtC9OnQ9u29v3o6AYjElWhexQaTR2if4dNjNKxlpyL3Pr1g5gY+wR1z54wYAC8+abfTbHZbBQXF5dxJqkOlfUo9KCuRqPR1AR3UVvHjoXmzSE1FQID7f70rVr5pXkRobi4mKKiIvLy8lyOIt27d/d5cqxqDT0ppWKUUn18aoFGo9HUFE/DPr7EYgGni3dKin04acMG91FbS0rsayC++MK+Hxdn92jyESUlJRQUFJCWlsb+/fvZv38/x44do6ioiLCwMK/czGtClUKhlPpKKRWtlGoB/AysVEot9Ys1Go1G4y2+SNZTWAiffgrHj9v3Dx+2uzs6h66//RZCQk6ttD1wAK67Dm6+2XPU1uJin61/sNlsGI1GTpw4wV9//cUff/zB4cOHOXnyJEFBQURFRREVFUVoaCgBAf6bcvam5mYikgdcAawUkbMBzz6XTYitW7fSqVOnCq6cNUUpxbRp01z7JSUltGrVivHjx1ernk6dOnHixIlal9FoGi2ekvV88QX88Yd96AegoADuuw+2bbPvp6TYJ5NXrbLvp6fD3/8On31m3xeBHTsgI8O+37UrPPzwqTUNQ4bYhWn1ar9EbRURzGYzOTk5HD16lN9//50///yTtLQ0bDYbERERREVFER4eXqcBP70RikClVBvgH8AWXzaulLpIKbVfKfWnUureSsoNUkpZlVJX+rL9yti6dSvjx4/nyJEjjB8/3idiERERwa+//orREQXys88+o127drWuV6M5ragsWc+ECXDmmfDqq/ZjQUHw9NP2xW0ALVrAuHH24SqAhATYvh0uvdS+36UL/PUXXHKJfb91a3jwQejWzb4fFgYdOtgXxvkoaqtzOCk1NZX9+/dz4MABkpOTXcNJUVFRREZGEhQUVG9ZEr0RikeAT4C/ROQHpVQX4GBtG1ZKGYAXgYuBs4DJSqmzPJR70mFDneAUCWdI7KKiIp+JxcUXX8wHH3wAwFtvvcXkyZNd506ePMnEiRPp06cPQ4cOZe/evQBkZWUxduxY+vfvzw033FBmHPLNN99k8ODB9OvXjxtuuOG0cIfVnL6ICHLttZUm65GWLbFcdhlWqxUJDgazGZzxy0JD4ZVXwBmBIDgYhg+HmsQIKx/i20uRsNlsFBUVVRhOys7Odg0nRUZG+n04qTpU6fUkIu8A75TaPwRM8kHbg4E/HfWhlFoHXAbsK1fuVuB/wCAftAnA/Pnz2bNnj9tz2dnZ/PrrrxUSChUVFTFmzBh69epFTExMhev69evHs1VFGwSuueYaHnnkEcaPH8/evXuZPXs227dvB2DRokX079+fjRs38uWXXzJ9+nT27NnDww8/zPDhw3nwwQf54IMPWL58OWB3tVy/fj3ffvstQUFB3HzzzaxZs4bp06dX7wvRaBogVqsVi8WCxWLBaDRSVFSE0WikQ6dORKak4O7Z2hYaypGnnqIwIAB+/x2lFEFBQYSEhJR5NxgMBAYGYjAYMBgMNX9Sd4rFrFn24SY3IiEiWCwWl3dSQUEBNpsNpRQhISFEREQ0+HzqVQqFUqo7sAyIF5FeDq+nCSLyr1q23Q44Vmo/GSgTHlEp1Q64HLiAKoRCKTUXmAvQoUOHGhu1f/9+j1nnbDYb+/fvZ+jQoTWuv0+fPiQlJfHWW29xibN76+Cbb77hf//7HwAXXHABWVlZ5Obmsm3bNjZs2ADAuHHjXEL1xRdfsHv3bgYNsn81RqORuLi4Gtum0dQHTv9/i8WC2WymsLAQo9F4Kpik1UqLTZtQ551HWMeO5Dz7LMadO4m9+24CSiXzsYWFkbFiBWrYMJw5EUUEq9XqilFmtVoRkQo35uDgYIKCgggODnaJiVNInO8eb+ZuoraWlJRgNBrJz88nLy/P9VkCAwMJCwtrMD0Fb/FmHcUrwELgZQAR2auUWgvUVijcfevlfbueBe4REWtViisiy4HlYF9wV1nZyp78yw87lSY8PJwtW7YwqpZZoyZMmMCCBQv46quvyMrKch1359rm/NzuPr+IMGPGDP7973/Xyh6Npi4QEUpKSiguLnbduI1GIyaTyXXzVkoRGBhIUFCQK6S9IS2N1v/+N3lZWeTccQfWVq0oGj+ejJYtiZszhwCj0SUSpmHDyrTprK8qu5xiZTKZsNlsrid+p10i4hKTkJAQgoODCQ4OLiMkxcXFFBYWkpub6/pMgYGBBAcHuz5LY8UboQgXke/L3ajcx42uHslA+1L7CUBKuTIDgXWOtmOBS5RSJSKy0Qftu2XUqFFs2bKlglj4SiQAZs+eTbNmzejduzdflUpwcf7557NmzRoeeOABvvrqK2JjY4mOjnYd/+c//8lHH33kyg89evRoLrvsMu644w7i4uI4efIk+fn5dGxAUSc1pydWq9XVSzCZTK5eQuneuvMm627oJfDwYcK3biVv9mysrVuTumkTxc4JZQemYcPIWLGC2IULObF4cQWR8BallGsIyhNOMSkpKcFisWC1Wl2fpbztjWU4qTp4IxQnlFJdcTztOzyPUn3Q9g9AN6VUZ+A4cA0wpXQBEens3FZKvQ5s8adIOCkvFr4UCYCEhARuv/32CscfeughZs2aRZ8+fQgPD2eVw4Vv0aJFTJ48mQEDBjBixAjX0NpZZ53Fv/71L8aOHYvNZiMoKIgXX3xRC4WmznCuDnbOJTgFwWw2l+kNO3sI3g65RG7cSPRrr1Fw6aXYWrWi2JFcqzymYcNI/uYbn30eT3gjJk2ZKmM9ObyclgPnANnAYWCqiLjJAVfNxpW6BPvwkgF4TUQeU0rdCCAiL5Ur+zp2oXi3qnp9Fetp69atzJo1i5UrV/pMJDSnN40p1pNzfL/0yzl05BymcQ7VOEXBOWxU7RtqSQlRa9di7t0bS//+qKIiAgoLsfop/EVTpaCggB49etQohEeNYz05XFNvEpExSqkIIEBE3GeeqQEi8iHwYbljL3koO9NX7XrLqFGjSPJTakGNpr6w2WxuBcDZK3D2EEpKSioMn4gIBoOBgIAADAYD4eHhPhliUWYzzV58kaJx4zjZvz8SHo7V04I2TZ1TqVA4JpHPdmwX1o1JGo2mujjH0MsLgPOmX/rdarWWubk7J2ydAhAQEOCagPXnOHvQwYNEvv022fffj0REkLppE9b4eL+1p6k53sxR/KSU2ox9LYVLLERkg9+s0mg0bhERcnNzKzz9FxcXuy3rvPE7RaAhLeIK2bOHqPXryZ88mZIuXbC2bl3fJmk84I1QtACysK9lcCKAFgqNpo45efIkycnJLpfMgIAAlwtmg/eyKSkheuVKStq3p+iiiyiYNImiMWOwuVnAqmlYeLMye1ZdGKLRaCqnuLiYtLQ0IiMjG6f3jVJEbNqEpVcvii66CAICtEg0ErxZmb2SigvhEJHZfrFIo9G4JS0tDaBRiUTQ/v00W7aMrCeeQEJDSVu7Fmkk6T81p/BmsHIL8IHj9QUQDRT406gGQ10kRdFovKCgoIDs7GzCG5knkCE7m7Bt2wg6cABAi0QjpUqhEJH/lXqtwR5uvJf/TatnfJEUpRwGg4F+/frRs2dP+vbty9KlSz3GlXKSlJTE2rVrq92W0WhkxIgRPo0mu2fPHj788MOqC5bjoYceYsmSJV6Xnzx5Mn369OGZZ56pdls1aa8yRo4ciXNNzuOPP+46npSURK9e7v8NFixYwJdffumT9sG+yvn48eN+90LyCSI0e/FFol97DQDT0KEkf/MNlj46MWZjpibuD92Amkfdawx4SopSS7EICwtjz549/Pbbb3z22Wd8+OGHPPzww5VeU1OheO2117jiiit8OkxRmVC4ArjVkrS0NL777jv27t3LHXfc4dU1vmq7KkoLRWXceuutPPHEEz5rNysrC4vFQnBwsM/qrC2hiYkkDB9OaGJi2RNKEbx3L8G//eY6JI2sF6SpiDepUPOVUnnOF/A+cI//TfMzI0fC66/bt4uL7ftvvmkXg3Hj3CdFufhi+/kTJ+zl33/ffs4xdlwd4uLiWL58OS+88AIiQlJSEueddx4DBgxgwIABfPfddwDce++9bN++nX79+vHMM894LFeeNWvWcNlll7n2n3rqKXr37k3fvn259157jqg9e/YwdOhQ+vTpw+WXX+6KITVy5EjuueceBg8eTPfu3dm+fTsWi4UHH3yQ9evX069fP9avX89DDz3E3LlzGTt2LNOnT+fIkSOMHj2aPn36MHr0aI4ePVrBrueff56zzjqLPn36cM0111Q4P3bsWDIyMujXrx/bt2+v1Mb777+fESNG8Nxzz1WoZ9++fYwcOZIuXbrw/PPPu457yt9x0003MXDgQHr27MmiRYsq1HfvvfdiNBrp168fU6dOBexP+tdffz09e/Zk7NixroRUHTt2JCsryzWnUBtMJhPp6elERETUui5fEZqYSNycOQQeP07cnDlEvP028dOnY8jMBCDz+ec58fTT9Wylxpd4M/QUJSLRpV7dReR/dWFcvTBrFpQKXVwGs9l+3kd06dIFm81GRkYGcXFxfPbZZ/z444+sX7+e2267DYAnnniC8847jz179riC/7krVxqLxcKhQ4fo5Mji9dFHH7Fx40Z27tzJzz//zN133w3A9OnTefLJJ9m7dy+9e/cu07spKSnh+++/59lnn+Xhhx8mODiYRx55hKuvvpo9e/Zw9dVXA7B79242bdrE2rVrmTdvHtOnT2fv3r1MnTrVrW1PPPEEP/30E3v37uWllyouwt+8eTNdu3Zlz549nHfeeZXamJOTw9dff81dd91VoZ4//viDTz75hO+//56HH36Y4uLiMvk79uzZg8FgYM2aNQA89thj7Nq1i7179/L111+7kkaVttvZI3Rec/DgQW655RZ+++03mjdv7goRDzBgwAC+/fZbN3917xERUlNTCQwMbDBrH5wi4QzvHWA00vLBBwnet4/Aw4fthUJC6tFCjT/wxuvpCxEZXdWxRkepqK0EBZ3ab9fOfZpFOJULNza27PW1WCjkjLVVXFzMvHnzXDewA47Jv/J4U+7EiRM0b97ctf/5558za9Ys10RoixYtyM3NJScnhxGOTF8zZszgqquucl1zxRVXAHD22WdXGsZkwoQJhIWFAZCYmOjKmzFt2jSXIJWmT58+TJ06lYkTJzJx4kSP9QJV2ugUK3eMGzeOkJAQQkJCiIuLIz09vdL8HW+//TbLly+npKSE1NRU9u3bR58qxtU7d+5Mv379gIrfU1xcHCkp5YMhV4/c3Fzy8/OJbiATwOVFwkmA2YxNKZTOrthk8SgUSqlQIByIVUrFcCp/RDTQtg5sqx+cGavKi0UNcuFWxaFDhzAYDMTFxfHwww8THx/Pzz//jM1m8xi//plnnqmyXFhYGCaTybXvLlFLVYQ4ngoNBkOlcwCVDYm4a/ODDz5g27ZtbN68mUcffZTffvutynwBNWk7pNRTrfMzeMrfcfjwYZYsWcIPP/xATEwMM2fOLPP9eduGsdQN1GQyuQS0JpSUlJCSktJgvJxCt28n7qabKoiEkwCTidiFC+skkqum7qmsP3sDsBv4G/CjY3s3sAl7ruumSw1z4VaHzMxMbrzxRubNm4dSitzcXNq0aUNAQACrV692jZ1HRUWRn38qDqOncqWJiYnBarW6bnZjx47ltddec+XXOHnyJM2aNSMmJsaVhnX16tWuJ3dPlLelPOeccw7r1q0D7HMkw4cPL3PeZrNx7NgxRo0axVNPPUVOTg4FBZ49rWtiY2WMHj2ad999l4yMDMD+PRw5coS8vDwiIiJo1qwZ6enpfPTRR26vDwoKchsqwx0HDhzw6BXlDenp6a7EN/VB0IEDxDz6KMrxGwo+eBAJCsLm4QHGFhbGicWL69JETR3iUShE5DlHPogFItK51KuviLxQhzbWD06x6NjRZyLhnAzt2bMnY8aMYezYsa6J05tvvplVq1YxdOhQDhw44Hpa7tOnD4GBgfTt25dnnnnGY7nyjB07lm8cT3cXXXQREyZMYODAgfTr18/lOrpq1SoWLlxInz592LNnDw8++GAVX8ko9u3b55rMLs/zzz/PypUr6dOnD6tXr64wyWy1Wrn22mvp3bs3/fv354477igzROaO6tpYGaXzd/Tp04cLL7yQ1NRU+vbtS//+/enZsyezZ8/m3HPPdXv93LlzXUNnlVFcXMyff/7JwIFuIzZXSVFREVlZWXXamwjIzKTZsmUEJicDEJiSQtSaNQQdPAhA3rXXcuzHH8l47TVs5XpKnrLLaZoOHvNRKKUuEJEvlVJXuDvfkIMC+iofRWPmp59+YunSpaxevbq+TTnteO+99/jxxx959NFHK5yr6ndos9n466+/XKk3/YUymwnfsoXi7t2x9O5N4NGjJIwYQebTT1N4xRVQXIyy2RA3E9Ol5yq0SDQs/JWPorKhJ2cf/1I3r/HVtkJTp/Tv359Ro0b5dMGdxjtKSkrcemJ5w8mTJzGZTH4RidDt2wndsQOwx+Rp+cADRDhcvEs6dODo99/bRQIgKMitSMCpFKQl7dppkThNqDLDXWPEU4/ib3/7W8Nf2appsogIf/zxh8cehcVi4cCBA4SHh/vEHTbo4EEMKSmYHPM6bf/+d6ytW5PuSLEbePQoJQkJ0EBcbzW1p84z3Cml7qysUhFZWm1L6pHQ0FCysrJo2bKlFgtNnSMiZGVlefRmc66ZcOaP8ERoYiKxCxdyYvHiCk/yASdOEPLbbxgdwtB86VJC9u61eyIpRcZ//4u1TRtX+ZIOTTvAgsZ3VOZSEeV47wEMAjY79i8FtvnTKH+QkJBAcnIymY7VoxpNXRMaGkpCQoLbc/n5+eTm5hIVFeX2PJSdG4ibM4eMZcsgOBjTkCEQEED0ypU0W76coz/+iERFkb1ggT18huPBqKRrV798Lk3Tp8qhJ6XUp8AkZ65spVQU8I6IXFQH9tUId0NPGk1DxWq1cuDAAQIDAz0OGbhb7CZBQajiYlI++ADLWWcReOwYAdnZWHr10sNJpyn1MZntpANgKbVvATpV2wqNRuOWjIwMrFZrtUQCQBUXI8HBBDjWhZS0b2+P0nqaiMTGjREMH55Aly4dGT48gY0bG048rKaGN6t5VgPfK6Xew+4scTnwhl+t0mhOE4xGI5mZmURGRnosE7twoccV0cpiIfaf/zztVkRv3BjB/fe3xGi0i+Lx44Hcf39LACZOLKxP05ok3qRCfUwp9THgXGY7S0R+8kXjSqmLgOcAA/CqiDxR7vxUTkWqLQBuEpGffdG2RlPfiAjHjx8nODi40gns7IULaXnvvQS4CSvSlFdEm82QlWXgxAkDmZn2d+dr/fpIl0g4MRoDWLIkRguFH/A2PsAeINVZXinVQUQqxpCuBkopA/ZQIBcCycAPSqnNIrKvVLHDwAgRyVZKXQwsB4bUpl2NpqGQnZ2N0WisdAKbkhKaP/ccxV27EnToUJmeRX0tdtu4MYIlS2JISTHQtq2VBQuyvb45V3bzP7UfQGamgbw897lUoqJsGI3uPRePHzeQnBxIQkLd5Cg5XfAmeuytwCIgHbBiDw4oQG1TVg0G/hSRQ4521gGXAS6hEJHSyRZ2AO5dRjSaRkZxcTGpqalVh+kIDCTr4YexxcQQkJ9f7yui3Q353HdfS3JyAujf3+zhxu/dzT821kpsrJXu3Ys55xwTrVpZXcdiY62OfRuhocLw4QkcP+7u9qU477wEBg40MWFCIePGFdKiReVZJJsCpcW7fXt4/HGoItJMtfDG6+lPYIiIZPmuWVBKXQlcJCLXOfanOdqZ56H8AuBvzvJuzs8F5gJ06NDh7CNHjvjSXI3Gpxw7dswVjNAdATk5BP3xB+ahQ13HNm6M4LvHfuapE9dxd+yrnPN/fWs1zCICJpOioEBRWBhAQUGAa9u+X3a7oCCAzZsjKgz5eKL0zf/Ujd7dvv3mXx3KCxZAWJiNBQuyMZkC2LQpggMHggkMFM4/38iECYVceGER4eFNb4Gxu+8iPByWL6+eWFTm9eSNUGwFLhQRn/bllFJXAX8vJxSDReRWN2VHAf8FhnsjWNo9VtOQKSgo4NChQ0RFRXlc/NnynnuI+PBDkrdvx9a8udubQUiIjZtuymXAALPjhq4cN/tT26WPl775FxQEUFSksFq9W3waHm4jIsJGZqaBUxkHSiO88kpGKRGwERLi35tyZUNgIvDHH0Fs2hTJ++9HkJISSHi4jQsvLOKyywoZPtxIDTxIGxxms2L48HacOFGxd9WxI1SSSqYCtRWKFdgX3X0AmJ3Ha7syWyk1DHhIRP7u2L/PUe+/y5XrA7wHXCwi7rP5lEMLhaahYrVa+fPPP1FKVRrPKSA3l+B9+1xDS56HWtwTGChERtqIiHC+27cjImxERtqIjDy1fep4xfMRETbCwwVn6nVPdrRrV8I33yRX78uoI2w2+OGHEDZtiuTDD8PJzTXQooWVceMKueyyQgYMMNNYgjXk5gawe3cIu3aF8MMPoezdG4LF4t54peyf3VtqKxQVEwgDIvKwu+PVMCoQOACMBo4DPwBTROS3UmU6AF8C08vNV1SKFgpNQyUjI4P09HT3E9gWC1Fr1pA/fTquO7ODzp074ulJ/u2308rd6IXgYPHLzc/TkM/jj2c1Cm8jiwW+/jqMTZsi+fzzMMzmABISipkwoZCJEwvp1s27fCN1RUqKgV27QvnhB7swHDgQhIgiMFDo3dvMwIFm/ve/SE6erDj348sehTfusQ87KokQEZ/9EkSkRCk1D/gEu3vsayLym1LqRsf5l4AHgZbAfx1d9BJPH0SjaeiYTCYyMjI8zkuEf/opLR95hOJu3TCVSvr0xhuevaLatbMyaJDZ43lf4xSDmno91TfBwXDhhUYuvNBIQYHik0/C2bw5kpdeasZ//9ucs84yM2FCIRMmFNKmTd1GXrbZ4M8/g/jhh1PCkJJiv0VHRNgYMMDMuHGFDBpkpm9fM2Fh9of8s86yuJ2jeOwx39nmTY9iGLACiBSRDkqpvsANInKz78zwLbpHoWloiAhJSUlVpkgN/uUXLL17A1BSAo880oLVq6M56ywzhw4FYTI1zif5hk5mZgBbtkSweXMke/aEoJQwZIjdc+qSS4po1sz3nlMWC/zySwg//BDCrl2h7N4dQk6OvWfQqlUJgwaZGTjQxODBZnr0sFBZssOKXk+q2l5PtR162glcCWwWkf6OY7+KSM3zPPoZLRSahkZOTg5Hjx4lOjq67AkRmj3/PIUTJ1LSsaPrcF6eYt68OLZvD+P663O5555s3n+/5usXNN6TlBTI5s0RbNwYyeHDQQQHCyNHFjFhQiGjRxur7aHlJC9P8dNPp3oLP/8cjNlsF/4uXYodomBi4EAzHTqU1GjosM7DjJdGRI6V887Q2XA0Gi8pKSkhJSXF7ZoJQ1oa0atWQUAAubfaHf6OHAnkuuviSEoK4oknTnD11fa84hMnFmphqAM6dSrhtttyufXWXH79NZhNmyJ4//0IPv00gqgoG3//u30SfNgwU6XinZ5ucInCrl0h/PFHMDabwmAQeva0cO21+QwaZObss03ExjbstR7e9CjeBZYCLwBDgduAgSJyjf/Nqxm6R6FpSBw/fpycnByPcxOG1FSsrVuDUnz/fQg33hiHCCxblsnQoRXDdmjqHqsVduwIZfPmCD76KIL8/ACioqwUFQWUcTEOChL69jWTnm7g2DH7U314uI3+/c0MGmRi0CAz/fqZ/bKeIzExkbvuuotVq1Zx4YUXVvv62g49xWKPxzQGe7TZT4Dbfb0Az5doodA0FIqKivjzzz8rrJmIWrkSgoPJLzWQ/O67kdx/f0vaty9hxYp0OnXSYSgaImaz4ssvw7jzztgyc0ZOAgKEsWOLXHMMZ55p8fuajcTERObMmYPRaCQ8PJwtW7YwatSoatVRK6FojGih0DQEbDYbf/31FyJSds2EzUbc9dcjwcFk/ve/2ETx1FMxvPxyM84918iLL2b6ZfJU41u6dOmISMWJBKWEQ4fqLjJEaZFwUhOxqFU+CqVUF6XU+0qpTKVUhlJqk1Kqi9etazSnKSdPnsRkMpUVCREICCBj2TIyn32WwqIAbrwxjpdfbsa11+axcmW6FolGQtu27qdqPR33B+5EAuw92fHjx7N161aftONN0Ja1wNtAG6At8A7wlk9a12iaKBaLhbS0tDLzEmFffEHcrFmowkIIDiYlK5x//KM1X3wRxqJFWTzyyMkmEVbidGHBgmzCwsqKujPeVF0gItx2220VRMJJUVERs2bN8klb3ng9KRFZXWr/TcdCOY1G4wYRITU1lYCAgDJ5JgLy8gjIywObjZ9/Dub66+MwGgNYsSKDkSPd/7NrGi71tfjwxIkT/O9//2PdunWcOHHCY7nw8HBWrlzpkza9mcx+AsgB1mEPL341EII9lwQictInlvgQPUehqU/y8vJISko6NYFtsdiXBANYrWz5KIoFC2Jp1crKihUZdO/esMJGaBoeNpuNb775hnXr1vH5559TXFzMwIEDueaaa2jZsiU333yzX+covOlRXO14v6Hc8dnYhUPPV2g0DqxWK8ePHycsLMwe+G/fPuLmziXzP//B1K8//3mxBc88E8PAgSZeeimDli31fITGM2lpabzzzju8/fbbJCcnExMTw/Tp07n66qvp1q2bq9yKFStq7fVUGd7Eeurss9Y0miZORkYGVqvVFabD2rw5xZ06UdQsjrvuiGXTpkiuuKKAxx8/QUhIPRuraZCUlJTw1VdfsW7dOrZu3YrNZuOcc87h7rvvZuzYsYS4+eEMGzaMFStWuNZR+FIkwLuhp0exhwO3OvajgedExDezJH5ADz1p6gOj0cjBgweJjIwkMD8fW3Q0KEVmpoG5c+PYsyeEhQuzuemm3EYT1lpTdyQnJ7N+/Xreffdd0tLSiI2N5aqrruIf//gHnTp18qoOf4Xw8MbrKRD4XinVRyk1Fns48N3VtkKjacKICMePHyc4OJjA3FzaXHopzZ99lt9/D2LixDbs3x/EsmUZ3Hxz7UQiMTGR4cOHk5iY6DvjNfVGcXExH330ETNmzOD888/nxRdfpEePHixbtozvvvuOu+++22uR8CfeDD3dp5T6AtgJZAPni8iffrdMo2lEZGdnYzQaiYqKwhYcTNG4cXwVeTHXXtWGyEgbb7+dRq9ellq1Udpnfs6cOaxYsYJhdZwzW+MbDh8+7Oo9ZGVl0aZNG2699VauuuoqEhIS6tu8ClQpFEqp87GH8HgE6A28oJSaLSIp/jZOo2kMFBcXk5qaSmRJCQFZWVhbtGRxy8d4/PEYevWy8MorGcTH124RVvmFVVosGh9ms5mPP/6YdevWsWPHDgwGAxdccAHXXHMNI0aMwGComHyooeCN19MS4CoR2QeglLoCe9a5v/nTMI2msZCWlgYitL71VlTWSWb1/o517zTn4osLefrpE64EMzXF0+pbLRYNg8TERBYuXMjixYvd/h0OHDjAunXreO+998jJyaF9+/YsWLCAK6+8kvj4+HqwuPp4M5ltcE5klzrWUgcF1Gjsk4eHDh0iKiqK4i928dpiK08euIZbbsnhzjtzCPBmFtADJSUl/PLLL8ycOZO8vDyP5dq1a8c333xT84Y0Naa0iIeFhblE22g08sEHH7Bu3Tp2795NUFAQY8eO5ZprruGcc84psxDTl9RnPoquSqllQLyI9FJK9QEmAP+qtiUaTRPCarVyPCmJ5vv2cbDVcK57bCIpKYEsXZrJ5ZdXf3WuzWbjjz/+IDExke+++47vv/+eggJ7LgqlFJ4e6oYMGUJqaipt2rSp1efRVA93w4GzZs3ivPPOY+fOneTn59O5c2fuu+8+Jk2aRMuWLevZ4prjTY/ia2Ah8LLOcKfRnCIjIwNZtIhWr7zK2WG/cizkDF5+OYOzz/Yuh7WIcOjQIZcw7Nixg+xse5ygzp07c8455zBs2DCGDh3KgQMHKgw/BQcHc+aZZ7J3716UUowePZqpU6dy3nnn+e2JVWPH03Cgk3PPPZdbb72VwYMHlwkv72/qs0cRLiLfl/uwOlC+5rTGZDKRkZHBh7F387PtPIztOrPx1VQSEir/10hOTnYJQ2JiIunp6QC0bduW0aNHM2zYMIYNG1ahd+BcUOVumOPYsWOsW7eO9evX89lnn9G+fXumTJnClVdeSWxsrN++g9OVtLQ0brnlFo8iAZCUlMSQIUPq0Cr/4k2P4iNgHvCOiAxQSl0JzBGRi+vCwJqgexQafyI2G+lLn+O+PVN4fU08o0YV8dxzmURFVfxfyszMLCMMR48eBaBly5auHsM555xDhw4dvHryrGzi1GKx8Omnn7J27VoSExMJCgrioosuYsqUKQwZMqROn2ybEikpKezYsYOdO3fy/fffk5SUVGn50iJe1/irR+GNUHQBlgPnYF9HcRiYKiJ1l5mjmmih0Pian5YupdXdd5P51FO07fA34q8ax2xWoGZP4v77s3F6Nubk5LBz506XMBw8eBCA6Ohohg4d6hKGbt26+fXG/ddff7F27Vreffdd8vLy6Nq1K1OnTuWKK66gWbNmfmu3KZCcnMzOnTtd4nDs2DHA/jccPHgwQ4YMYejQoeTk5DB37twyPYv6FAmoR6EoVUkEECAi+dW2wHOdF2Ffo2EAXhWRJ8qdV47zlwBFwEwR+bGqerVQaHzJq1OXMnntXUQAhcA1EU9iNg5g5CN9mXBZOrt27XIJw2+//YaIEBYWxuDBg13CcNZZZ/nET15EXJPalf3vOs+ZTCY+/PBD3nrrLfbs2UNISAjjxo1j8uTJ9OnTp4JYVXU/8HS+Npkya3qt87qgoKAa3RiddRw7dswlCjt27CAlxb5ErHnz5gwePJihQ4cyZMgQevToUeFv6Mnrqb6od6HwNUopA3AAuBBIxh4aZLJzvYajzCXArdiFYgj2GFNVDvxpodD4itIi4aQQuClhJPvi8/n5558pKSkhODiYAQMGuIShT58+ZTPblUJEsNls2Gw213bpd0+ICEop181KKeW2V1L6mHNbKcW+fftYt24dmzdvprCwkJ49ezJ58mQmTJjgSrDk7lpvt93hycbKbK7suLtyubm5WCwWIiIiqpzEFxGSkpJcovD999+TmpoKQIsWLcr0GLp37+6VU0BV6yjqkqYoFMOwBxv8u2P/PgAR+XepMi8DX4nIW479/cBIEUmtrG4tFBpf8NPSpXS/q6xIOCkEbmjblqgJExg6dCh9+/YlJCTEJQBw6qbm/B8rfdMODAzEYDBgMBgIDAx0vZzHDAaDK/FR+Vdtyc/PZ82aNSxbtoy9e/cSFRXFtddey4033kifPn1qXX9dsnXrVmbOnMkzzzxDt27dMBgMrhDvcMqzbOfOnS5xyMjIAOzzRM7ewpAhQ/w+HFgXNEWhuBK4SESuc+xPA4aIyLxSZbYAT4jIN479L4B7RKSCCiil5gJzATp06HD2kSM1m0Kx2WzatVADQHJgIAlWz6E3jhkMmP/4g6CgINdN3nnD93STNxgMDeZmJCLs3LmTl156ifXr12MymRg2bBg33XQTV155pStUekNl69atjB8/nqKiIsLDw9mwYQM9evTgp59+Yu/evfzwww/s3LnTlQUuLi7OJQpDhw6lS5cuDeZv4Svq0z0WpdQ5QKfS5UXkjWpbUq5aN8fKq5Y3ZZz2LMc+6c7AgQNrpH4Wi4Xk5GQ6dOhAYKBXX42miVBSUoLFYsFsNlNUVERhYSF/3nknsYsXE+qmfCFw4qmn6H/GGXVtqs9QSjF06FCGDh3K0qVLWbVqFS+99BLTp09n/vz5zJw5kxtuuIHu3btXuHbr1q3MmjWLlStX+jz3gTeUFgmw54ceN24ckZGR5ObmAnZhOPfcc13i0Llz5yYnDHWFN0EBVwNdgT2A8/FKgNoKRTLQvtR+AlA+0KA3ZXyGiJCbm0tKSgrt27fXP6omiIiUEYXCwkIKCwspLranI3UOCwUGBvJX+4UE8D5j+aNMPP5C4K0pT3PdnXfWy2fwBy1atOCOO+5g/vz5fPXVV7z00ks8//zzLF26lAsuuIAbb7yRyy67jODg4DI36fHjx/ssm5rVaiUnJ4eTJ0+6XtnZ2RX2Dxw4wPfff+8a4it9fUFBAQsWLOD666+nefPmZGRkEBAQUGY4SlN9vHGP/R04S3w8RqWUCsQ+mT0aOI59MnuKiPxWqsw47Gs4nJPZz4vI4Krqrukchdls5uDBg9hsNhISEmjRokW169A0HESE4uJiLBYLJpPJJQpWx3CSUxSCgoIqeLNs2RLObbe1JCpsCH8v3sfKYqPL6+mtKU9z3ZqmIxKeSEtLY+XKlbz88sscOXKE+Ph4xowZw4YNGyrNz2w0Giu90Xvad/YEPBEVFUWLFi1ISUlxCbs7Onbs6FrrYDabSUtLIy8vj9DQ0Bp7RzUW6nMdxTvAbVVNINcEh1fTs9jdY18TkceUUjcCiMhLDvfYF4CLsLvHznI3P1Ge2gpFWFgYRUVFdOvWjdBQdwMPmoaGiGCxWLBYLBiNRgoLCzEajS5RCAgIcIlCVXNQ21ZnEbroKf7Vfhi/HL2DJ554gjEGA/H33kvmU0/Rvwn1JLzBarXyySef8K9//ctjwiSlFDExMRQVFWEymTzWZTAYaNGiBS1atCAmJsa1XdV+8+bNXTe/8sNOpXGXL1pEKCws5Pjx4157RzVW6lwolFLvYx9iigL6Ad8DriA2IjKh2pbUEbUVisjISMxmMwEBAXTp0qVBx4k/HbHZbC5RcM4nGI1Gl/uoUso1wVzdG8JHH4XzybwdrAiYw6XhxRi7dObLL78kJibGT5+m8dCpUycqcxKJjIzkpptuqvTmHxkZ6ZMhIHdi4U4kSmOz2cjKyiI9Pb3JDkfVx2T2kmq31IQICQmhoKCAtLQ02rVrV9/mnNaYzWaMRiNFRUUUFRWVGfYwGAwEBQURERFR63/6jz8K47bbWtG3/1juOeMidry9hg2LFtG8efNafoKmwcqVKyt9kt+8eXOdTWyPGjWKLVu2lPF6qmquJCAggFatWtGsWTNSU1NPm+EoX+DxcUtEvhaRr4FLnNulj9WdifVHREQEJ06cqHLsVOM/8vPz+fPPPzl27Jjr7xAZGUlUVBRRUVGEh4cTFBRUa5H4Ygt0umUmczp9zP33b+P1d9/iyiuvZMyYMU3uqbOmOG/O4eHhZY57c5P2pz0dO3asVvvBwcF06NCBzp07Y7PZyM/PrzAxrimLN/3yC90ca7ABAX2JUoqIiAiOHTuG2exd6GiN78jNzSUpKYng4GCioqIICwsjMDDQ5zfuTz8N48H5oXQOSeHeGw6zePEiIiMjuffee10rljV2yotFfYlEaXuSkpKq3b5SisjISLp160br1q1dvdX6WlfW0PEoFEqpm5RSvwA9lFJ7S70OA3vrzsT6xbmQKjk5WT911CHZ2dkcOXLE1WPwF59/GsK8W1rRqlcM5m83sCXUxo4dO5g3bx5nnnmm7k24oaZP8g0R53BU9+7diYiIID8/v1KPqtOVyiazmwExwL+Be0udyheRk3VgW43xxWR2efLz84mLi2s0OW4bKyJCVlYWKSkpRERE+NWR4MsvQrHNfZDo5vC3L/+PwKBCxowZQ/Pmzfnwww/p1KmT39rWNEwKCgoatXeUvyazK/sWDEAecAuQX+qFUuq0W2AQERFBeno6+fk+C56rKYeIkJGRQUpKCpGRkX4Via++CuOmm+KgZTNGTAwiupnw4osvkpqayr333kvbtm391ram4eIcjmrTpk2jGI6y2WyYzWYKCgooKCjwWw+4Mq+n3ZwKl1G+dQG6+MWiBorTne7YsWN069ZNe0r4GBEhNTWVEydOEBkZ6dcnua+/DuP2uVF071HMoDdvxRhtJSkpiVdffZVLL72UsWPHeoz8qmn6BAQEEBsbS3R0NGlpaeTk5BAaGlrvvwmbzUZxcXGZobGAgAAiIiJo2bIloaGhhISE+OXe5FEoRKSzz1tr5AQFBVFcXMzx48fp2LGjHr/2ETabjePHj5OdnU1UVJRfv9dt20L5ds4Wfgl4nPyn3yKiWRygePTRRwkKCuLOO+/U6UM1wCnvqBYtWnD8+HHy8/PrbDjKW1Hwh3OHO7wNChgDdINT8dFEZJu/jGrIhIeHk5eXR1ZWlr6h+ACr1cqxY8fIz8/3u0hs3x7K3LlxTGzfg8gze2Pr1ByAL774gi+//JIFCxbQu3dvHRBSUwbncNTJkydJS0vz+WK98qIgIhgMhnoTBXd4ExTwOuB27AH59gBDgUTgAr9a1oCJjIwkJSWF8PDwCj7lGu8pKSnh6NGjFBUVERUV5de2vv02lHuvC6BzlxLuXtOBvBbPAnYHhkcffZQuXbowbdo0vQJb4xZfDUc1BlFwhzePTrcDg4AdIjJKKfU34GH/mtWwCQgIIDQ0lGPHjtG1a1f9BFoDiouLSUpKori42K2XmS/57rtQ/jP7OPtLRpEx60kCWpxaBvTqq69y5MgRXn75ZTp27NjovFw0dUv54aiCggLCw8Pd/m4aqyi4w5s7nElETI4YOiEi8odSqoffLWvgBAcHU1hYSGpqKgkJCQ3+D92QMJvNJCUlYbPZ/N4j27EjlDlz4ujRMRRjn8sIHDUQ52qYlJQUXnzxRcaOHcsFF1xAdHS0X23RNB3KD0c5IxE3BVFwhzdCkayUag5sBD5TSmXjx5wQjYnw8HCys7OJiIjQIcm9xGQycfjwYZRSfs+gtmNHCA/MstClXSGvrs2jMPbRMucff/xxbDYb8+fPp02bNo3yH1hTf5QejsrIyMBqtTYJUXBHlUIhIpc7Nh9SSm0FmgEf+9WqRoIzxMfx48cJDw/XIcmroKioiMOHDxMUFOR3V8OdO0O4eVY0e0p6E9G1N/mxz5U5n5iYyAcffMC8efM466yzdKgOTY0JDg4mISGhvs3wK956PQ0HuonISqVUK6AdcNivljUSnNFLjx07pkOSV0J+fj5JSUl1Eq3zhx9CmD07njbtSiiZOg/zgK5lzhcXF/PQQw+RkJDAtGnT9Gp7jaYKvPF6WgQMBHoAK4Eg4E3gXP+a1ngIDQ2loKCA9PR0vaLXDbm5ua64Tf6e+N+1K4RbZ4Zzbotf+dfaGALirsBSrsybb77JgQMHeO6552jXrp3uCWo0VeCNi8flwATsWSARkRTsyYw0pdAhyd2TlZXFkSNHiIiI8LtI7N4dwsyZ8SxTN/NB4QXEh+VUKJOZmckzzzzDeeedx6hRo4iLi/OrTRpNU8Cb/1yLiIhSSgCUUnow1w1KKcLDw0lOTiYsLKzel/vXNyJCZmYmaWlpfg/JAfDTT3aRaNXKSqelt5Cdei7iZm3G4sWLMZlMLFy4kPj4eB2KRaPxAm/+e99WSr0MNFdKXQ98DrziX7MaJ87Um6d7SHIRIS0trQ5FIpiZ01txVfj7vLU2lZj+bSm6pGJurZ9++ol33nmHmTNn0rVrV1q2bOlXuzSapkKV/8EisgR4F/gf9nmKB0XkP/42rLESFhZGYWEhmZmZ9W1KvWCz2UhJSSEzM5OoqCi/i8SePcHMmNGamSFreC1jIp3+/NqjXQ899BBxcXHMnj2bNm3aaMcDjcZLvBo0FpHPgM/8bEuTwRmSPCIiwu+rjhsSNpuN5ORkcnNz/R63CWDv3mCmT29N8+ZWrl47koy9L2IaPtxt2bfffpu9e/eyZMkSYmJiaNasmV9t02iaElU+7imlrlBKHVRK5Sql8pRS+UqpvLowrrFSOiT56ZIty2q1cuTIEfLy8oiMjPS7SPzySzDTpsUzJWg9by/7nbYJYh9uctNubm4uixcvZuDAgVx44YW0adNGh+rQaKqBN/8tTwETRKSZiESLSJSI1CrWgVKqhVLqM4cAfeaITlu+THul1Fal1O9Kqd+UUrfXps26JigoyDUM05ATn/iC4uJiDh8+TFFRUZ2IxK+/2kWia3gKLxbO5swNz1VafunSpeTk5PB///d/REZG+j0AoUbT1PBGKNJF5Hcft3sv8IWIdAO+oGyqVSclwF0icib2iLW3KKXO8rEdfiUiIoLc3FyysrLq2xS/YbFYSEpKwmw218nq5t9+C+baa+OJiBCee1uRvu4tcu6+22P533//nTfffJMpU6bQpUsXWrdu3WTCKmg0dYXHOQql1BWOzV1KqfXYYz2ZnedFZEMt2r0MGOnYXgV8BdxTuoCIpAKpju18pdTv2FeE76tFu3VOZGQkqampTTIkudls5vDhw4hInYjEvn1BXHttPBcEbuPh25OJan8elvb9PJYXERYtWkSzZs246aabiImJaXJ/A42mLqhsMvvSUttFwNhS+wLURijiHUKAiKQqpSpd9aSU6gT0B3ZWUmYuMBegQ4cOtTDNtwQEBBASEtLkQpIbjUYOHz7smo/xBxs3RrBkSQwpKQZatbJSUBBAs2grK9r8k4jXskiZ9AFU4rn0/vvv88MPP/Cvf/2LqKgovbhOo6khlaVCnVWbipVSnwOt3Zz6v2rWE4ndNXe+iHicRBeR5cBygIEDBzaoSYHg4GAKCgqaTEjywsJCDh8+THBwsN8WFm7cGMH997dkiPFrvmEWszJW8hUjueOOPPL+sYyCoqJKRaKwsJDHH3+c3r17c+mllxIbG0tISIhfbNVomjp+e7wVkTGeziml0pVSbRy9iTZAhodyQdhFYk0th7rqnYiICLKzs4mMjGzUWdTy8vI4evSo35K4O1myJIYhxq/ZwngiKOIDxrGZCdyz8nWuuy4aaxW5I1544QXS09N54YUXXOGgNRpNzagvH8HNwAzH9gxgU/kCyv7YvQL4XUSW1qFtfqF0SHKTyVTf5tSInJycOosA2/34NpdIAIRj5GrWc37K21Vee+jQIVasWMGkSZPo0aMHrVu31qE6NJpaUF9C8QRwoVLqIHChYx+lVFul1IeOMucC04ALlFJ7HK+KcRkaEQaDgcDAQI4dO4bVaq1vc6rFiRMnOHr0aJ0E9wtNTGQLl7pEwokCXlY3EZqY6PFaEeGRRx4hNDSUu+66i8DAwEbdg9NoGgLehBkPASYBnUqXF5FHatqoiGQBo90cTwEucWx/g/3e0KRobCHJ6zq4H0DU7XcTXk4knIRLEcELF5L8zTduz3/++ed8/fXX/POf/yQyMlKH6tBofIA3//WbsLuzlmAPNe58aWqIMyR5Xl7DXuBe18H9wJ6Z7qr81ylS7t1YbWFhnFi82O05k8nEo48+Srdu3bjmmmsIDQ3VoTo0Gh/gzRhCgohc5HdLTiOcIcmPHTtGt27dGlRIcqvVislkorCwkJycHEwmU53EbQL49ttQrr8+jrbtWnLkztfosWAWAUaj67wtLIyMFSswDRvm9vrly5dz7Ngx3nzzTaxWq86DrdH4CG8eEb9TSvX2uyWnGQ0lJLnNZqOoqIgTJ05w6NAhfv/9dw4dOkRmZiYBAQFER0fXyc3266/DmDMnjg4dSnhrbSotIkxkvPoqNscajapEIjk5mWXLlnHxxRdz9tlnExUVpfNgazQ+whuhGA7sVkrtV0rtVUr9opTa62/DTgecIclPnDhRZ22KCCaTiezsbI4cOcLvv//OX3/9RVpaGiUlJURERLhusnW1OPDzz8OYOzeOrl2LWbs2jQ57PiF+5kwCCgrIWLGCknbtKhUJgMceewyA+++/H4vFokN1aDQ+xJs7wcV+t+I0JiIigrS0NCIiIvz2BGyxWDCZTOTl5ZGXl+fyuAoODiYsLKxeI6l+9FE4t93WirPOsrBqVTrNm9swjhlD5vPPUzRmDAQEeJy4dvLtt9/y8ccfc+eddxITE0Pz5s39tlpcozkdqSzWU7RjJXR+Hdpz2uEMgXH06FHOOOMMn/j7l5SUYDKZKCgoIDc3F4vFglIKg8FAaGhogwmxvXlzBHfeGUu/fmZeey2d6EgrqqAIiYyk8NJLq64Ae+Tahx56iA4dOnDddddhtVpp1aqVny3XaE4vKutRrAXGA7uxx3Yq3Y8XoIsf7TqtCAoKwmKxkJKSQocOHao9ZGKz2VwT0Lm5uZhMJkQEg8FAcHBwgwyr/b//RXD33bEMHGhmxYp0IiOFqJWraPbqq6Ru2IA1Pt6relatWsWff/7JK6+84hKJhuQcoNE0BSqL9TTe8d657sw5fSkdkryqcBMigtlspqioiNzcXAoLCxERAgICCA4OJiIiokGPz69bF8n997fknHNMLF+eQXi4PTSXpV8/isaMwepl8L7MzEyee+45Ro4cyciRI7FYLDoPtkbjB5pGKNMmgjMkeURERJkxdhHBYrFgNBrJz88nPz8fq9WKUqpRCENp3ngjikWLWjJiRBEvvZRJaOip+I3m/v0x9+/vdV1PPPEEZrOZBx54AJPJRJs2bZpMdF6NpiGh/6scrFkD998fzLFjPWnb1sqCBdlMnFi36wqdIcmPHj1Kp06dsFgs5Ofnk5eX50qpGhQU1KDmGarDq69G89hjLRgzpogXXsjAGcw1+pVXUCYTubfcAl5+rt27d7NhwwZuuukm2rdvj81m06E6NBo/oYUCu0jMnQtFRfan8uPHA7n/fvsQRl2LRXBwMIWFhRw8eBDANc8QGhpap3b4mv/+txmLF8dw8cWFPPdcJq45exGC9+9HFRW5zXftDqvVykMPPUTr1q255ZZbMBqNdOzYsVGKp0bTGPBKKJRSBiCesrGejvrLqLrm//4PisqFFjIaA1iyJKbOhQJoUgvFROD555vx7LMxXHZZAUuWnKDM6JBSnFiyBMxmr4Vi/fr1/Prrrzz//PMEBgYSEhJCdBVhxzUaTc2p8hFMKXUrkA58BnzgeG3xs111ylEPkpeSooPJ1QYRWLKkOc8+G8OVV+bz9NNlRSJi0yYMaWn2HS+TCmVnZ7N48WKGDBnCuHHjMJvNOlSHRuNnvOmr3w70EJGeItLb8erjb8PqEk+ZU9u0aVyhwBsSIvDYYzH897/NmTw5nyefzCqTkC4gN5eWDzxAsxdeqFa9S5cuJS8vj0WLFmE2m2nWrFmT6oFpNA0Rb4TiGJDrb0Pqk8ceg/AKwUqFyEgbRqN+Uq0uNhssWtSCFSuaMWNGHo89llVhjtrWrBkpmzaRc/fdXte7b98+1q5dy7Rp0/jb3/5GSUkJ8V6ut9BoNDXHG6E4BHyllLpPKXWn8+Vvw+qSqVNh+XLo0EFQSmjXroTJk/M5eDCI666Lc01ya6rGZoP772/J6tXRXH99LosWnaww9RB4+DAAJZ07Y/NybkFEWLRoEc2bN+eOO+7AaDTSokWLRj/Jr9E0BrwRiqPY5yeCgahSrybF1Klw4ICFvXt/45tvknn88ZM8/fQJduwIZfbseAoLtVhUhdUKd98dy/r1UdxySw733ZddQSRCfvyRdmPGEPH++17VmZiYyPDhw1m8eDG7du3i7rvvJioqCpvNRpyXC/M0Gk3tqNLrSUQeBlBKRdl3pcDvVjUQLr+8EIMB7rgjlpkz41m50h5qQlORkhK4665YNm+O5I47srntNvejlZazziJn/nyKLrigyjoTExOZM2cORqORZcuW0aVLF6666iqKioqIi4vTebA1mjrCG6+nXkqpn4Bfgd+UUruVUj39b1rDYMIEu9//Tz+FMGNGPPn5umdRHosFbr21FZs3R3LPPSc9igQ2GxIaSu6ttyJVTECXFgknx48f57vvvsNgMOhQHRpNHeLN0NNy4E4R6SgiHYG7gFf8a1bDYvz4Iv7zn0z27g1h+vTW5OVpsXBiNsPNN8fx8ccRPPDASW680X1619Bt22gzcSKGlJQq63QnEva2zFx//fUcPHhQ58HWaOoQb4QiQkS2OndE5CvgtPNHvPjiIl58MYPffgtm2rTW5ObqVcAmk+KGG+L44otwHn00i9mzPecAVzYbtogIbF70BBYuXFhBJE61aeLWW2+tsc0ajab6eOX1pJR6QCnVyfH6J3DY34Y1RMaONbJsWQZ//BHMtdfGk5Nz+opFUZFizpw4tm0L44knTnDttZWnLTGOHEn62rWIFwvrFi9e7DFUeFhYGK+//npNTNZoNDXEmzvdbKAVsAF4z7E9qzaNKqVaKKU+U0oddLx7jOamlDIopX5SSjWI1eCjRxt56aUMDhwIZurUeE6ePP3EoqBAMXNmPDt2hLJkyQmuvtqzf0PYF18Q+e679hV4Xq6ezszMdEXHLU1oaChbtmxh1KhRtbJfo9FUjyrvciKSLSK3icgAEekvIreLSHYt270X+EJEugFfOPY9cTvwey3b8ymjRhlZvjydv/4KYsqU1mRlnT5ikZenmDEjnh9/DOG55zK54orKY2FFbthA1MqVdreoKhARli9fzu23387AgQNZvny5K9x6aGgoGzZs4AIvvKU0Go1v8XiHU0o963h/Xym1ufyrlu1eBqxybK8CJnqwIQEYB7xay/Z8zogRJl59NYMjRwKZMqU1mZlNXyxycgKYNq01v/wSwosvZjJ+fFGV12T+5z+kv/46VOHKarVaeeSRR/j3v//NuHHjWLVqFWPGjGHFihW0bduW119/nYsv1unbNZr6oLJ1FKsd70v80G68iKQCiEiqUsrTyqlngbtpoAv8hg838dprGcyZE8fkya1ZuzaduLimGR/q5MkApk2L588/g1m2LIPRo91PNjsJ+eknLGecgURFYasih7XZbGb+/Pl8/PHHzJkzh/vvv98VMvzss8/m008/5YwzzvDZZ9FoNNXD42OwiOx2bPYTka9Lv4B+VVWslPpcKfWrm9dl3himlBoPZJSyo6ryc5VSu5RSuzIzM725xCcMG2bitdfSSU0N5JprWpOW1vTcNjMzA5g8uTV//RXEK6+kVykSymgk7rrriL3vvirrzsnJYdq0aXz88cf885//5J///KdLJCwWCxaLhS5duhDiZXRZjUbje5RI5SuNlVI/isiAcsd+EhHvc1ZWrHM/MNLRm2gDfCUiPcqV+TcwDSgBQoFoYIOIXFtV/QMHDpRdu3ZV2y6z2czBgweJjIys9rU//BDCrFnxtGplZe3atEYdeXbjxgiWLIkhJcVAfLwVqxUKCgJ49dUMzjnH5FUdIT/+iDUmhpLOnlOuHz9+nJkzZ3L06FGefvppxo8f7zpXXFyMxWKhc+fOhFeM2KjRaHyMUmq3iAx0d66yOYrJSqn3gc7l5ie2Alm1tGkzMMOxPQPYVL6AiNwnIgki0gm4BvjSG5GoLwYNMvPGG+mcOGHgmmtac/x44+xZbNwYwf33t+T48UBEFGlpgWRmGpgzJ9crkVAFdg8o84ABlYrEvn37mDRpEunp6axataqCSJjNZi0SGk0DobIZ2O+Ap4E/HO/O113ARbVs9wngQqXUQeBCxz5KqbZKqQ9rWXe9MWCAmTfeSCM72y4WycmNL9PskiUxGI3lfxaK996repooMDmZhPPPJ2Jz5b4O3377LVdffTVKKd555x2GDh3qOldcXIzJZNIiodE0ICqbozgiIl+JyLBycxQ/ikjVvo6VICJZIjJaRLo53k86jqeIyCVuyn8lIuMr1tTw6N/fwptvppGXF8DVV7fm6NHGIxZFRcpjT8ibbH+2iAiKRo/GPGCAxzIbN25k1qxZtGvXjg0bNtCjx6kRR6dIdOnSRScj0mgaEN4EBRyqlPpBKVWglLIopaxKKc+xGjT06WNhzZp0iooUkye3JimpYYtFaqqBJ5+M4ZxzEgD3i+Latq16zsUWE0PW4sWUJCRUOCcivPTSS9xxxx2cffbZvP3227Rp08Z1vqSkBLPZTKdOnbRIaDQNDG+c/18AJgMHgTDgOuA//jSqKdCrl4U1a9IwGu1icehQwxOLn34K5rbbYjnvvASWL4/mnHNM3HZbNmFhtjLlwsJsLFjgeY1l4KFDtLrhBgwZGW7PW61WHnroIZ588knGjx/P66+/TnSphEUlJSUYjUY6duxIVFSD9ITWaE5rvLp7icifSimDiFiBlUqp7/xsV5PgrLOKWbs2jWuvbe1aZ9G1a3G92lRSAh9/HM5rr0Xz00+hREXZmD07j+nT80hIsPcaOncucXk9tW1rZcGCbCZO9LwCO3j/fkJ++cUepqMcJpOJ+fPn88knn3Dddddx3333udxf7fZokdBoGjreuMduA8ZgXx2dBqQCM0Wkr//Nqxn14R5bGQcP2kN9AKxdm0a3bnUvFrm5Abz1ViSrV0eTkhJIp07FzJyZx6RJBb5JxmQ2Q7m1Djk5OVx//fXs3r2bBx54gFmzyoYIs1qtFBYW0qlTpzI9DI1GU/fUyD22FNMAAzAPKATaA5N8Z17Tp1u3YtatS0MpYfLk1uzfX3eZ2f76K5AHHmjBsGEJPPlkCzp1KubVV9P54ovjzJiRXyuRCNq/n9Bt2+w75UQiOTmZK6+8kr179/Kf//zHo0h07NhRi4RG08DxJhXqEcemEXjYv+Y0Xbp2Leatt9KYOrU1U6a0ZvXqNM46yz89CxH45ptQVq6MZuvWcIKDhYkTC5g5M48zz/Rdm81feIHQHTtI/vprpJQr6759+5g5cyZms5k33niDIUOGlLmutEg0a9bMZ/ZoNBr/4FEolFJvi8g/lFK/ABUeO0Wkj18ta4J07VrCunVpTJnSmqlTW/Pmm+n07GnxWf0mk+K99yJ4/fVoDhwIJjbWyh13ZDNlSj6xsbaqK6gmJ5YsIfDQoTIisX37dm6++Waio6N588036d69e5lrnCLRoUMHLRIaTSOhsh7F7Y73RrF+obHQqZNdLCZPbs2UKfGsXp1Onz61E4v0dAOrV0exdm0U2dkGzjrLzJIlmYwfX1h+RMgnGI4fxxoXh4SEUHzmma7jGzZs4J577qFr1668/vrrtG7dusx1TpFo3749zZs3971hGo3GL1S24C7VsXkFUOJYgOd61Y15TZMOHUpYvz6N6Ggb117bmj173Gdzq4q9e4OZPz+W4cMT+O9/mzFokIl161LZsiWVSZP8IxIUFxM/cyZxN9/sOiQiLFu2jLvuuotBgwbx9ttvVxAJm81GYWEhCQkJxMR4zFOl0WgaIN64x0YDnyqlTgLrgHdFJN2/ZjV9EhJKXHMW06e35vXX0xkwwFzldSUl8NlndvfWXbtCiYy0MX16HjNm5NOhQ60WzHtHUBA5CxZgc7iyOtdIvPnmm0yYMIGnnnqqQqRXm81GQUEB7dq1o0WLFv63UaPR+JQq3WNdBZXqA1yN3eMpWUTG+NOw2tDQ3GMrIyXFwNSprcnMNPD66+kMHOheLPLyAli/PpJVq6I5fjyQ9u2LmTkzn6uuyicqygfurd5gs0GpNRAmk4nbb7+dTz/9lLlz53LPPfeUWSNhv+SUSLRs2bJu7NRoNNWmMvfY6iwXzsC+jiIL8JRoqFGjvvqK7rNmcfLppzENG1YnbbZta+Wtt+wT3DNmxDN7di7vvRflWuw2Y0YuyclBvPtuJEVFAQwZYuLBB08yenQRhjoMUKsKC2l99dXkzptH0UUXkZ2dzfXXX8+PP/7Igw8+WMH9FU6JRNu2bbVIaDSNGG8W3N2EvSfRCngXWC8i++rAthpTox7F1q3I+PGooiJsYWFkrFhRZ2IBkJFh4NJL25CRYaBsvCUhIAAmTixk1qw8evXynZdUdTBkZBA7fz458+fzV9u2zJgxg+TkZJ599lm3KUpLi0RsbGw9WKzRaKpDbXsUHYH5IrLHp1Y1JLZuBYdIAAQYjcTNmVOnYhEXZ3WM6pQPyqeIiyvh6adP1IkdnrDGxZG+Zg2/7dvHrCuuwGKx8OabbzJo0KAKZZ0i0aZNGy0SGk0ToMqV2SJyLxCplJoFoJRqpZTynJGmseEQCRwi4cQpFqGJiXVmSnq6+7EkT8f9SWhiIgnDhxP2+efEPPYYKj+fbdu3c/XVVxMUFMS7777rUSTy8/Np3bq1FgmNpolQZY9CKbUIGAj0AFYCQcCbwLn+Na2OmDWrgkg4CTAaiV24kORvvqkTU9q2tXL8eMU/iTchvn1JaGIisTNnEmixEHvTTQQAmyMjmfPCC3Tr1o2VK1cSHx9f4brSPYlWrVqhlPuQ5RqNpnHhTayny4EJ2OM8ISIpQNMJ87lyJXjIpGYLCeHEU0/VmSkLFlQ/xLevKS0SAIaSEiwivPbsswwZMoT169e7FQkRoaCggPj4eC0SGk0TwxuhsIh9xlsAlFJNK6vMqFGwZUsFsbCFhKCsVoL37q0zUyZOLOTxx7No164EpYR27Up4/PGsSkN8+5LyIuEkxGrlw4AA3po7120ocBEhPz+f+Ph44uLitEhoNE0Mb7yeFgDdsOe2/jcwG1grIg02eZFPvJ5eeYWgAwcovPRSbLGx9pVugQ0v+ZAvaTVoEBEnPE+aF8bGkvnDD2WOOUWiVatWtG7dWouERtNIqVWYcRFZgt0t9n/Y5ykebMgiUWNGjaJ4wwYsbdrYvZ3OPZf8WbPsIiFC3M030+Lhph0893mLBfezNfZxx/IrJbRIaDSnB95muPsM+MzPttQ7MnIkBz/9tOLKbKuV4s6dsboZm28qZO3cyd15ebwB/AMoPb5YCEwKDmbG88+7jjnnJGJjY7VIaDRNnMrCjOfjJry4ExE5fbLNBAaSfd99rt2QHTsI++47cm69FYLqLgmRr1EmE9nvv8+/d+1iw4YNjDYYKDn3XNYlJrKhuJgISonE668zzLGmxNmTaNmyJW3atNEiodE0cTwKhYhEASilHsEeumM19tVgU2lKXk81IGzbNiI++ojcG29EGqlQ/PHHHxTNm8f4v/7ih+Bgpk6dyty5c2nbti2JiYlMmjmTlywWbnQjEgUFBbRs2ZK2bdtqkdBoTgO8mczeKSJDqjpWrUaVagGsBzoBScA/RKSCD6hSqjn2XN29sPduZotIlSvg6iIoYEBeHrboaCgpIfzjjykaNw4awU3zwEcfsWrNGtZ++y0dw8O57YILOO/BB2nVqlWZcomJiSxcuJDFixdXEImYmBjatWunRUKjaULUNoSHVSk1FXuIcQEmA7VdAXYv8IWIPKGUutexf4+bcs8BH4vIlUqpYMD9god6wObI8xzx/vu0uvNO0po3xzR8eD1b5R4RYefOnbzy/POsS0ykMCiIuPnzmTFjhscEQsOGDeObUgsNnSLRvHlz3ZPQaE4zvBGKKdhv2M9hF4pvHcdqw2XASMf2KuArygmFUioaOB+YCSAiFqB+IuJVQuHEiVhbtHCJhCE9vcFMeosIX3/9NVufeoo3fv+d2NhYPrjqKobeeCP9u3SpVl0FBQU0a9aMdu3aVQglrtFomjZVCoWIJGG/sfuSeGcGPRFJVUq5C1veBcgEViql+gK7gdtFxO3qM6XUXGAuQIcOHXxsbiUohWnECMAuEm3//nfyrr+e3FtuqTsbymGz2fj000958cUX6f3rr6wFRs2YwdB77yU0NLTa9eXn59OsWTMSEhK0SGg0pyF+W0GmlPocaO3m1P95WUUgMAC4VUR2KqWewz5E9YC7wiKyHFgO9jmK6ltce6wtWpA3ezaFzrDbInU6b1FSUsL777/Pay++SOFffyGdOnH2v/5FenExI6dOrZGHVkFBAdHR0VokNJrTGL8JRWUZ8JRS6UqpNo7eRBvsSZHKk4w9k95Ox/672IXCr9hsNkSkZmPwQUHk3naba7fFQw8hoaFk33uvXwXDbDazYcMGXnrpJY4ePcrOsDA6t21L7iefEBgcjLEGdToD/DVr1oz27dtrkdBoTmPqKybFZmAG8ITjfVP5AiKSppQ6ppTqISL7gdGAXxMmBQcH06pVK06cOEFQUFCNhmlc2Gz2l9XqN5EwGo2sW7eO5cuXE5KWRkzv3vzfyy/TPjQUa3AwgcHBNa63pKSENm3a0LJlSy0SGs1pjtdCoZQaCjwOhACLRWRjLdp9AnhbKTUHOApc5WijLfCqiFziKHcrsMbh8XSIilEkfIpSijZt2tC8eXOOHz9OXl4eERERGGqSczQggJOPPmoffgKCDh4kbOtW8ubMobY5TPPz81m9ejWvvfYaWVlZTO7Vi9VZWWRPnkzB2LG4z7pdNcXFxZhMJqKjo2ndujUhISG1slOj0TQNKluZ3VpE0koduhN7uHEFfAdsrGmjIpKFvYdQ/ngKcEmp/T3Yc2HUKWFhYXTt2pWcnBxSU1MREcLDw2s2HOW4JuK994hav56CSZOw1TB/dHZ2NitXruT1118nPz+fy4cOZfKyZQwaOJC8//wH4wUX1Khem81GUVERBoOBDh06EB0drd1fNRqNi8rGFF5SSj2glHKOv+Rgd4u9Gsjzt2H1jVKKmJgYunfvTvPmzcnPz8dsrumzOuQsXEjK5s12kRCxZ87zsNgxMTGR4cOHk+jIrpeRkcFjjz3G8OHD+c9//sM555zDgUmTeOfPPxncowcoRe5tt9XILddkMrliNnXv3p1mzZppkdBoNGWoLITHRKXUpcAWpdQqYD52oQgHJtaJdQ2AwMBA2rVrR0xMDCkpKeTn5xMeHl794SilsLZrB0DYF18Qf/31ZLz8MkVjx5YplpiYyJw5czAajcyaNYvzzjuPbdu2UVJSwqRLLmHu3Lmc0bs3wXv3ktujBxIWVqPPZbVaKSwsJDIykk6dOtVuPkaj0TRpvAnhYQBuBsYBj4nI9rowrDbUNIRHVdhsNrKzs0lLs4/I1Xg4ymolYtMmCi+7DAwGVH4+EhVVRiRKM3LkSB5ZuJChN95I4aWXkrNwYY0/g4hQVFSEUoq2bdvqHoRGowFqmI9CKTVBKfUN8CXwK3ANcLlS6i2lVFf/mNqwCQgIoGXLlmWGoyyWGiwWNxgovOIKl0i0HTeODEdIDaPRyEjgMPal6xHAzp07Sc7NpWDSJEznnVdj+81mM/n5+TRv3tz1GbRIaDSaqvDYo1BK7QWGAWHAhyIy2HG8G/CoiFxTZ1ZWE3/1KMpTWFhISkoKJpOpRsNRhYWFfPnhh7R+/nmeSU4mEbs4bMEuECagGOgDFLdrVyb2UnWwWq0UFRURGhpKu3btCPeQI1yj0Zy+1DQoYC72XkQYpRbEichBx/HTnoiICLp27Up2djapqakopaocjrJYLGzbto3Nmzfz2WefYTKZaNu2LWdfeiljP/6Y94uLca5+CAUMwNCgICYsXlxt+5zDTCJC27ZtiYmJ0WsiNBpNtalMKC7HHim2mNoHAWyyOIejoqKiSE9PJzs7m9DQUIJLLXaz2Wx8//33bN68mQ8//JDc3FxiYmK48sormTBhAmeffTbhO3cS+8knFf4gQcBqpTiBvYfhLcXFxRiNRmJiYoiPjy9jj0aj0VSHKiezGyN1NfTkjoKCAlJSUjCbzSQlJfH+++/z/vvvk5aWRnh4OBdeeCGXXXYZw4cPJ6hU7KWE4cMJPH7cY70l7dqR7MXQk81mo7CwkODgYNq1a+dVbg2NRqOpbT4KTTVIS0tj3bp1vPnmmxw8eJDAwEBGjBjB/fffz+jRoz3OD5xYvJi4OXMIMFaMzGQLC+OEF0NPztAb8fHxxMbG6mEmjUbjE7RQ+IC0tDTWr1/P2rVr+f777wEYMWIE8+fP59xzz0UpRVhYWJkeRHlMw4aRsWJFBbGwhYWRsWIFJkeWOXc4h5mio6Np06aNDr2h0Wh8ihaKGpKbm8uGDRtYu3YtX375JTabjf79+7N48WKuvvpq2rdvD9gnlAsLCzl+/Dj5+flERER4fNIvLxZViUTp0BsdO3bUoTc0Go1f0HMU1cBkMvHBBx+wdu1aPvjgA8xmM127dmXKlClMnjyZM8880+O1NpuNrKws0tLSCAwMJDQ01ONNPTQxkdiFCzmxeLFHkTCZTFgsFlq1akVcXFzNAhdqNBqNg8rmKLRQlGLr1q3MmjWLlStXMmrUKMCeDGjr1q2sXbuWDRs2kJeXR3x8PNdccw1Tpkxh0KBB1XqKN5vNpKWlkZubW+VwlDtKh95o06YNYTUM4aHRaDSl0ZPZXrB161bGjx9PUVER48ePZ/Hixezfv5/169eTnp5OdHQ0kyZNYsqUKYwcOZLAwJp9dSEhIXTo0IH8/HyXd1R4eHiVE8/ONREA7du316uqNRpNnaGFgrIiAVBUVMQtt9xCUFAQl156KVOmTGHcuHE+C5ynlCI6OpqIiAiysrJIT08nMDDQY+/AYrFgMplo2bIlcXFx1e6FaDQaTW047YWivEiUJjAwkHnz5rmGoXyNwWAgLi6OZs2akZKSQl5eXpnhqNKhN8444wwdekOj0dQLp/0cRadOnThy5IjH8x07diQpKclHlnlGRMjLyyMlJYWSkhICAgIQEVq3bk2LFi30mgiNRuNXahQ99nRh5cqVHp/Uw8PDWblyZZ3YoZSiWbNmdO/e3dXL6N69u144p9Fo6p3T/g40atQotmzZUkEswsPD2bJli9+GnTxhMBiIj48nISFBx2fSaDQNgtNeKKCiWNSXSGg0Gk1DRAuFA6dYdOzYUYuERqPRlOK093oqzahRo+pk4lqj0WgaE/XSo1BKtVBKfaaUOuh4j/FQ7g6l1G9KqV8dKVh9s5BBo9FoNF5TX0NP9wJfiEg34AvHfhmUUu2A24CBItILe7I3nVlPo9Fo6pj6EorLgFWO7VXARA/lAoEwpVQgEA6k+N80jUaj0ZSmvoQiXkRSARzvceULiMhxYAlwFEgFckXkU08VKqXmKqV2KaV2ZWZm+slsjUajOf3wm1AopT53zC2Uf13m5fUx2HsenYG2QIRS6lpP5UVkuYgMFJGBrVq18s2H0Gg0Go3/vJ5EZIync0qpdKVUGxFJVUq1ATLcFBsDHBaRTMc1G4BzgDeranv37t0nlFKe43I0DmKBE/VtRANBfxdl0d9HWfT3cYrafBcdPZ2oL/fYzcAM4AnH+yY3ZY4CQ5VS4YARGA14FcBJRBp9l0IptctT3JXTDf1dlEV/H2XR38cp/PVd1NccxRPAhUqpg8CFjn2UUm2VUh8CiMhO4F3gR+AXh63L68dcjUajOX2plx6FiGRh7yGUP54CXFJqfxGwqA5N02g0Gk05dAiPhovuPZ1Cfxdl0d9HWfT3cQq/fBdNMh+FRqPRaHyH7lFoNBqNplK0UGg0Go2mUrRQNCCUUu2VUluVUr87giHeXt821TdKKYNS6iel1Jb6tqW+UUo1V0q9q5T6w/EbGVbfNtUnp3vQUKXUa0qpDKXUr6WOeRVwtbpooWhYlAB3iciZwFDgFqXUWfVsU31zO/B7fRvRQHgO+FhE/gb05TT+XnTQUABeBy4qd6zKgKs1QQtFA0JEUkXkR8d2PvYbQbv6tar+UEolAOOAV+vblvpGKRUNnA+sABARi4jk1KtR9c9pHTRURLYBJ8sd9jbgarXQQtFAUUp1AvoDO+vZlPrkWeBuwFbPdjQEugCZwErHUNyrSqmI+jaqvqhu0NDTiCoDrtYELRQNEKVUJPA/YL6I5NW3PfWBUmo8kCEiu+vblgZCIDAAWCYi/YFCfDSs0BipbtBQTe3QQtHAUEoFYReJNSKyob7tqUfOBSYopZKAdcAFSqkqA0I2YZKBZEdoG7CHtxlQj/bUN66goSJSDDiDhp7upDsCrVJJwNVqo4WiAaGUUtjHoH8XkaX1bU99IiL3iUiCiHTCPkn5pYictk+MIpIGHFNK9XAcGg3sq0eT6htX0FDH/81oTuPJ/VI4A66C54Cr1aa+osdq3HMuMA34RSm1x3HsfhH5sP5M0jQgbgXWKKWCgUPArHq2p94QkZ1KKWfQ0BLgJ06zUB5KqbeAkUCsUioZe1y8J4C3lVJzsIvpVT5pS4fw0Gg0Gk1l6KEnjUaj0VSKFgqNRqPRVIoWCo1Go9FUihYKjUaj0VSKFgqNRqPRVIoWCk2doZSyKqX2OKJ9vq+Ual6DOgYqpZ73cC5JKRVbQ9seUkotcHP8RqXU9JrUeTqg7HzpiEVVWbnPfRXJVFP3aKHQ1CVGEenniPZ5EriluhWIyC4Ruc33pnls7yUReaOu2vMHjqB5/uIS4GcvQs2sBm72ox0aP6KFQlNfJOKIjKuU6qqU+lgptVsptV0p9TfH8ascvY+flVLbHMdGOnNTKKVaKqU+dQTJexlQjuOdysXoX6CUesixfb1S6gdHnf9TSoVXZmTpnoZS6iul1JNKqe+VUgeUUuc5jhuUUkuUUr8opfYqpW51HB/tsO0XR+6AEMfxJKXU40qpRKXULqXUAKXUJ0qpv5RSN5Zqe6HD1r1KqYc92DfHYctXSqlXlFIvOI6/rpRaqpTaCjyplOqnlNrhqOs959O947qBju1YR8gUlFIzlVKbHH+X/UqpRR6+oqmUWv2rlLrW8f3sUUq9rJQyOE5tBiZX9l1rGi5aKDR1juPmMRr7zQPsK2pvFZGzgQXAfx3HHwT+LiJ9gQluqloEfOMIkrcZ6OBF8xtEZJCjzt+BOdU0P1BEBgPzHe0DzMUenK6/iPTBvno6FHu+gKtFpDf2KAg3larnmIgMA7Y7yl2JPQfJIwBKqbFAN2Aw0A84Wyl1fmlDlFJtgQcc110I/K2crd2BMSJyF/AGcI/Dvl9K2V4Zg7ELQT/gKqeglONcYLfDnjOBq4FzRaQfYHVcj4hkAyFKqZZetKtpYGih0NQlYY7QJFlAC+AzZY+Uew7wjuPcy0AbR/lvgdeVUtdjT0xTnvOBNwFE5AMg2wsbejl6Lb9gv4n1rOZncAZq3A10cmyPAV4SkRKHLSeBHtiD1h1wlFnlsNeJUyR/AXaKSL6IZAImx9zNWMfrJ+xhKv6GXThKMxj4WkROOgLjvVPu/DsiYlVKNQOai8jXHmzxxGcikiUiRsfnHu6mTAtH7hSwi//ZwA+Ov+Vo7OHRnWRgj/SqaWToWE+ausQoIv0cN64t2OcoXgdyHE+gZRCRG5VSQ7AnL9qjlKpQBnAXg6aEsg9BpVNkvg5MFJGflVIzscfKqQ5mx7uVU/8/yo0dyst6bKW2nfuBjuv/LSIvV1JHVW0UVnEeyn5X5VOJlv9Mbr9rpVSAiNgc9qwSkfs8tBUKGL2wSdPA0D0KTZ0jIrnY01guwH7jOKyUugpcXjR9HdtdRWSniDwInADal6tqG46hDaXUxYDTqyYdiHPMYYQA40tdEwWkKns496k++kifAjc6J42VUi2AP4BOSqkzHGWmAV97uN4dnwCzHT0ulFLtlFLlk9B8D4xQSsU42p7kriLH953tnFMpZ0sS9l4A2Ie/SnOhsudgDsOeKe1bN9Xv51Sv4QvgSqedjms7OrYV0NrRnqaRoYVCUy+IyE/Az9hDiE8F5iilfgZ+w56QBmCxYyL4V+yi8HO5ah4GzldK/Yh9mOaoo+5i7GP9O7H3XP4odc0DjuOflTteG151tL3X8RmmiIgJe3TXdxzDXDbgJW8rdGRrWwskOq5/F7vIlS5zHHgc++f5HHvY8VwPVc7A/n3uxT7n8Ijj+BLgJqXUd0B51+JvsHsr7QH+JyK73NT7AY5emYjsA/4JfOpo5zNODSOeDexwDs9pGhc6eqxG04hRSkWKSIGjR/Ee8JqIvOeDemcCA0VkXhXl2gBviMiFVZR7DtgsIl/U1jZN3aN7FBpN4+Yhx8Txr8BhYGNdNu7Iy/yKqmLBHfCrFonGi+5RaDQajaZSdI9Co9FoNJWihUKj0Wg0laKFQqPRaDSVooVCo9FoNJWihUKj0Wg0lfL/1+Ewv/c3uF8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(table.index,table['mean'],'-o',color='b',label='Data (no controls)')\n",
    "ax.fill_between(table.index, table['low'], table['up'], color='grey', alpha=.25)\n",
    "ax.plot(table.index,table['sim'],'-D',color='k',label='Model')\n",
    "ax.plot(table.index,table['mean_h'],'-D',color='r',linestyle=':',label='Data (controls for health)')\n",
    "plt.xlabel('Residual income group (e)')\n",
    "plt.ylabel('% deviation in health expenditures')\n",
    "plt.legend()\n",
    "plt.savefig('../figures/fig_d2_spending_by_income.eps', bbox_inches='tight',dpi=1200) \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Production Function Investigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 488,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wmean(x,var,wvar):\n",
    "    xx = x.loc[~x[var].isna(),:]\n",
    "    names = {var: (xx[var] * xx[wvar]).sum()/xx[wvar].sum()}\n",
    "    return pd.Series(names, index=[var])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 499,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                 totinc   R-squared:                       0.036\n",
      "Model:                            WLS   Adj. R-squared:                  0.035\n",
      "Method:                 Least Squares   F-statistic:                     31.77\n",
      "Date:                Fri, 27 May 2022   Prob (F-statistic):          9.34e-172\n",
      "Time:                        15:38:11   Log-Likelihood:            -2.8739e+05\n",
      "No. Observations:               24585   AIC:                         5.748e+05\n",
      "Df Residuals:                   24555   BIC:                         5.751e+05\n",
      "Df Model:                          29                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "===================================================================================\n",
      "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept        4.047e+04    850.982     47.559      0.000    3.88e+04    4.21e+04\n",
      "C(age)[T.51.0]    926.9075   1030.140      0.900      0.368   -1092.229    2946.044\n",
      "C(age)[T.52.0]   1947.4740   1041.384      1.870      0.061     -93.702    3988.650\n",
      "C(age)[T.53.0]   2113.8869   1045.247      2.022      0.043      65.139    4162.634\n",
      "C(age)[T.54.0]   2911.0148   1049.184      2.775      0.006     854.550    4967.479\n",
      "C(age)[T.55.0]   2024.0943   1069.727      1.892      0.058     -72.636    4120.825\n",
      "C(age)[T.56.0]   2273.5078   1101.792      2.063      0.039     113.929    4433.087\n",
      "C(age)[T.57.0]    972.2476   1111.962      0.874      0.382   -1207.264    3151.760\n",
      "C(age)[T.58.0]    681.0707   1140.463      0.597      0.550   -1554.305    2916.446\n",
      "C(age)[T.59.0]    -56.4336   1166.415     -0.048      0.961   -2342.679    2229.811\n",
      "C(age)[T.60.0]   -650.2533   1181.808     -0.550      0.582   -2966.668    1666.161\n",
      "C(age)[T.61.0]  -2780.4536   1214.208     -2.290      0.022   -5160.375    -400.533\n",
      "C(age)[T.62.0]  -3067.7348   1215.603     -2.524      0.012   -5450.391    -685.079\n",
      "C(age)[T.63.0]  -5347.4268   1286.441     -4.157      0.000   -7868.929   -2825.925\n",
      "C(age)[T.64.0]  -4967.5074   1297.281     -3.829      0.000   -7510.257   -2424.758\n",
      "C(age)[T.65.0]  -5552.5456   1259.615     -4.408      0.000   -8021.467   -3083.624\n",
      "C(age)[T.66.0]  -7086.4103   1293.591     -5.478      0.000   -9621.928   -4550.893\n",
      "C(age)[T.67.0]  -9296.4477   1344.193     -6.916      0.000   -1.19e+04   -6661.747\n",
      "C(age)[T.68.0]  -1.022e+04   1336.897     -7.647      0.000   -1.28e+04   -7603.007\n",
      "C(age)[T.69.0]  -1.264e+04   1361.326     -9.286      0.000   -1.53e+04   -9973.409\n",
      "C(age)[T.70.0]  -1.233e+04   1396.368     -8.829      0.000   -1.51e+04   -9591.147\n",
      "C(age)[T.71.0]  -1.293e+04   1398.046     -9.251      0.000   -1.57e+04   -1.02e+04\n",
      "C(age)[T.72.0]  -1.317e+04   1360.634     -9.680      0.000   -1.58e+04   -1.05e+04\n",
      "C(age)[T.73.0]  -1.314e+04   1377.066     -9.542      0.000   -1.58e+04   -1.04e+04\n",
      "C(age)[T.74.0]  -1.254e+04   1452.861     -8.630      0.000   -1.54e+04   -9689.821\n",
      "C(age)[T.75.0]  -1.332e+04   1451.900     -9.175      0.000   -1.62e+04   -1.05e+04\n",
      "C(yr)[T.2001.0]   381.2067    624.495      0.610      0.542    -842.840    1605.254\n",
      "C(yr)[T.2002.0]   636.1179    607.146      1.048      0.295    -553.924    1826.160\n",
      "C(yr)[T.2003.0]   784.1494    628.737      1.247      0.212    -448.213    2016.511\n",
      "C(yr)[T.2004.0]  1246.7679    626.265      1.991      0.047      19.250    2474.286\n",
      "==============================================================================\n",
      "Omnibus:                    13467.506   Durbin-Watson:                   1.487\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):           139437.632\n",
      "Skew:                           2.452   Prob(JB):                         0.00\n",
      "Kurtosis:                      13.586   Cond. No.                         23.1\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      m   R-squared:                       0.019\n",
      "Model:                            WLS   Adj. R-squared:                  0.018\n",
      "Method:                 Least Squares   F-statistic:                     16.60\n",
      "Date:                Fri, 27 May 2022   Prob (F-statistic):           1.51e-82\n",
      "Time:                        15:38:11   Log-Likelihood:            -2.6287e+05\n",
      "No. Observations:               24585   AIC:                         5.258e+05\n",
      "Df Residuals:                   24555   BIC:                         5.260e+05\n",
      "Df Model:                          29                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "===================================================================================\n",
      "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept        2650.8104    313.879      8.445      0.000    2035.589    3266.032\n",
      "C(age)[T.51.0]     92.3453    379.960      0.243      0.808    -652.399     837.090\n",
      "C(age)[T.52.0]     33.3876    384.107      0.087      0.931    -719.486     786.261\n",
      "C(age)[T.53.0]    378.9225    385.532      0.983      0.326    -376.744    1134.589\n",
      "C(age)[T.54.0]    300.3913    386.984      0.776      0.438    -458.121    1058.904\n",
      "C(age)[T.55.0]    659.2502    394.562      1.671      0.095    -114.114    1432.615\n",
      "C(age)[T.56.0]    761.2672    406.388      1.873      0.061     -35.279    1557.813\n",
      "C(age)[T.57.0]   1390.5152    410.139      3.390      0.001     586.617    2194.413\n",
      "C(age)[T.58.0]   1809.7532    420.652      4.302      0.000     985.250    2634.256\n",
      "C(age)[T.59.0]   1553.5968    430.224      3.611      0.000     710.331    2396.863\n",
      "C(age)[T.60.0]    834.5356    435.902      1.915      0.056     -19.858    1688.929\n",
      "C(age)[T.61.0]   1772.9062    447.852      3.959      0.000     895.089    2650.724\n",
      "C(age)[T.62.0]   1606.6983    448.367      3.583      0.000     727.872    2485.525\n",
      "C(age)[T.63.0]   2400.0795    474.495      5.058      0.000    1470.041    3330.118\n",
      "C(age)[T.64.0]   2427.8759    478.493      5.074      0.000    1490.000    3365.752\n",
      "C(age)[T.65.0]   3006.3361    464.600      6.471      0.000    2095.691    3916.981\n",
      "C(age)[T.66.0]   2583.3226    477.132      5.414      0.000    1648.115    3518.531\n",
      "C(age)[T.67.0]   2470.2066    495.797      4.982      0.000    1498.415    3441.998\n",
      "C(age)[T.68.0]   2642.3964    493.105      5.359      0.000    1675.880    3608.913\n",
      "C(age)[T.69.0]   2768.6491    502.116      5.514      0.000    1784.472    3752.826\n",
      "C(age)[T.70.0]   3962.1003    515.041      7.693      0.000    2952.589    4971.611\n",
      "C(age)[T.71.0]   3873.5205    515.660      7.512      0.000    2862.796    4884.245\n",
      "C(age)[T.72.0]   4648.8220    501.861      9.263      0.000    3665.145    5632.499\n",
      "C(age)[T.73.0]   4285.8670    507.921      8.438      0.000    3290.310    5281.423\n",
      "C(age)[T.74.0]   4205.2848    535.878      7.847      0.000    3154.932    5255.638\n",
      "C(age)[T.75.0]   5115.5602    535.523      9.552      0.000    4065.902    6165.218\n",
      "C(yr)[T.2001.0]   403.9923    230.341      1.754      0.079     -47.489     855.474\n",
      "C(yr)[T.2002.0]   747.6726    223.941      3.339      0.001     308.734    1186.611\n",
      "C(yr)[T.2003.0]  1023.6433    231.905      4.414      0.000     569.095    1478.192\n",
      "C(yr)[T.2004.0]  1152.6648    230.994      4.990      0.000     699.903    1605.426\n",
      "==============================================================================\n",
      "Omnibus:                    35343.528   Durbin-Watson:                   1.801\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):         18565481.207\n",
      "Skew:                           8.477   Prob(JB):                         0.00\n",
      "Kurtosis:                     136.553   Cond. No.                         23.1\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                  adl1p   R-squared:                       0.006\n",
      "Model:                            WLS   Adj. R-squared:                  0.004\n",
      "Method:                 Least Squares   F-statistic:                     4.778\n",
      "Date:                Fri, 27 May 2022   Prob (F-statistic):           3.51e-16\n",
      "Time:                        15:38:12   Log-Likelihood:                 16103.\n",
      "No. Observations:               24585   AIC:                        -3.215e+04\n",
      "Df Residuals:                   24555   BIC:                        -3.190e+04\n",
      "Df Model:                          29                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "===================================================================================\n",
      "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
      "-----------------------------------------------------------------------------------\n",
      "Intercept           0.0078      0.004      2.100      0.036       0.001       0.015\n",
      "C(age)[T.51.0]     -0.0007      0.004     -0.152      0.879      -0.009       0.008\n",
      "C(age)[T.52.0]      0.0018      0.005      0.396      0.692      -0.007       0.011\n",
      "C(age)[T.53.0]     -0.0029      0.005     -0.636      0.525      -0.012       0.006\n",
      "C(age)[T.54.0]   7.641e-05      0.005      0.017      0.987      -0.009       0.009\n",
      "C(age)[T.55.0]      0.0005      0.005      0.114      0.909      -0.009       0.010\n",
      "C(age)[T.56.0]      0.0016      0.005      0.339      0.735      -0.008       0.011\n",
      "C(age)[T.57.0]      0.0074      0.005      1.534      0.125      -0.002       0.017\n",
      "C(age)[T.58.0]      0.0062      0.005      1.256      0.209      -0.003       0.016\n",
      "C(age)[T.59.0]      0.0054      0.005      1.054      0.292      -0.005       0.015\n",
      "C(age)[T.60.0]      0.0017      0.005      0.335      0.737      -0.008       0.012\n",
      "C(age)[T.61.0]      0.0061      0.005      1.159      0.246      -0.004       0.016\n",
      "C(age)[T.62.0]      0.0017      0.005      0.321      0.748      -0.009       0.012\n",
      "C(age)[T.63.0]      0.0052      0.006      0.924      0.356      -0.006       0.016\n",
      "C(age)[T.64.0]      0.0055      0.006      0.969      0.333      -0.006       0.017\n",
      "C(age)[T.65.0]      0.0080      0.005      1.459      0.145      -0.003       0.019\n",
      "C(age)[T.66.0]      0.0110      0.006      1.952      0.051    -4.4e-05       0.022\n",
      "C(age)[T.67.0]      0.0147      0.006      2.508      0.012       0.003       0.026\n",
      "C(age)[T.68.0]      0.0070      0.006      1.202      0.229      -0.004       0.018\n",
      "C(age)[T.69.0]      0.0111      0.006      1.877      0.061      -0.000       0.023\n",
      "C(age)[T.70.0]      0.0159      0.006      2.612      0.009       0.004       0.028\n",
      "C(age)[T.71.0]      0.0210      0.006      3.457      0.001       0.009       0.033\n",
      "C(age)[T.72.0]      0.0233      0.006      3.927      0.000       0.012       0.035\n",
      "C(age)[T.73.0]      0.0313      0.006      5.214      0.000       0.020       0.043\n",
      "C(age)[T.74.0]      0.0360      0.006      5.699      0.000       0.024       0.048\n",
      "C(age)[T.75.0]      0.0243      0.006      3.842      0.000       0.012       0.037\n",
      "C(yr)[T.2001.0]    -0.0028      0.003     -1.030      0.303      -0.008       0.003\n",
      "C(yr)[T.2002.0]     0.0048      0.003      1.798      0.072      -0.000       0.010\n",
      "C(yr)[T.2003.0]     0.0059      0.003      2.143      0.032       0.001       0.011\n",
      "C(yr)[T.2004.0]     0.0006      0.003      0.215      0.830      -0.005       0.006\n",
      "==============================================================================\n",
      "Omnibus:                    31865.848   Durbin-Watson:                   1.732\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):          3487254.406\n",
      "Skew:                           7.614   Prob(JB):                         0.00\n",
      "Kurtosis:                      59.324   Cond. No.                         23.1\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "df = pd.read_stata('../data_sources/meps/MEPS0004_agghealth.dta', convert_categoricals=False)\n",
    "df = df[(df.age>=50)&(df.age<=75)]\n",
    "df = df[df.totinc>10e3]\n",
    "df = df[~df.mepsexp.isna()]\n",
    "data = df[['totinc','adl1p','age','mepsexp','totmcd','yr','perwt']].dropna(axis=0)\n",
    "data['m'] = data['mepsexp']\n",
    "nqs = 10\n",
    "for v in ['totinc','m','adl1p']:\n",
    "\tmodel = wls(v+' ~ C(age) + C(yr)',data=data).fit()\n",
    "\tprint(model.summary())\n",
    "\tdata[v] = data[v].mean()+model.resid\n",
    "data['qy10'] = pd.qcut(data['totinc'],q=nqs,labels=[x for x in range(1,nqs+1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 500,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_inc = data.groupby('qy10').apply(wmean,var='totinc',wvar='perwt')\n",
    "table_m = data.groupby('qy10').apply(wmean,var='m',wvar='perwt')\n",
    "table_h = data.groupby('qy10').apply(wmean,var='adl1p',wvar='perwt')\n",
    "table = pd.concat([table_m,table_inc,table_h],axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 501,
   "metadata": {},
   "outputs": [],
   "source": [
    "table['lhs'] = -np.log(table['totinc']-table['m'])\n",
    "table['rhs'] = np.log(table['adl1p'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 502,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>m</th>\n",
       "      <th>totinc</th>\n",
       "      <th>adl1p</th>\n",
       "      <th>lhs</th>\n",
       "      <th>rhs</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>qy10</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5704.409489</td>\n",
       "      <td>9736.314013</td>\n",
       "      <td>0.030254</td>\n",
       "      <td>-8.301994</td>\n",
       "      <td>-3.498117</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5062.608382</td>\n",
       "      <td>15612.693063</td>\n",
       "      <td>0.019061</td>\n",
       "      <td>-9.263889</td>\n",
       "      <td>-3.960131</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5425.270781</td>\n",
       "      <td>20183.170152</td>\n",
       "      <td>0.023251</td>\n",
       "      <td>-9.599534</td>\n",
       "      <td>-3.761391</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5363.768712</td>\n",
       "      <td>23866.689340</td>\n",
       "      <td>0.018706</td>\n",
       "      <td>-9.825684</td>\n",
       "      <td>-3.978934</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>4862.176647</td>\n",
       "      <td>27926.048720</td>\n",
       "      <td>0.013739</td>\n",
       "      <td>-10.046023</td>\n",
       "      <td>-4.287484</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>4654.088022</td>\n",
       "      <td>32897.991195</td>\n",
       "      <td>0.011841</td>\n",
       "      <td>-10.248633</td>\n",
       "      <td>-4.436195</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>5059.624238</td>\n",
       "      <td>39191.088204</td>\n",
       "      <td>0.008488</td>\n",
       "      <td>-10.437975</td>\n",
       "      <td>-4.769057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>4872.726798</td>\n",
       "      <td>47733.198064</td>\n",
       "      <td>0.009946</td>\n",
       "      <td>-10.665705</td>\n",
       "      <td>-4.610572</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>5016.744444</td>\n",
       "      <td>60975.390602</td>\n",
       "      <td>0.004505</td>\n",
       "      <td>-10.932368</td>\n",
       "      <td>-5.402655</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>4672.357166</td>\n",
       "      <td>104956.863777</td>\n",
       "      <td>0.005643</td>\n",
       "      <td>-11.515766</td>\n",
       "      <td>-5.177292</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                m         totinc     adl1p        lhs       rhs\n",
       "qy10                                                           \n",
       "1     5704.409489    9736.314013  0.030254  -8.301994 -3.498117\n",
       "2     5062.608382   15612.693063  0.019061  -9.263889 -3.960131\n",
       "3     5425.270781   20183.170152  0.023251  -9.599534 -3.761391\n",
       "4     5363.768712   23866.689340  0.018706  -9.825684 -3.978934\n",
       "5     4862.176647   27926.048720  0.013739 -10.046023 -4.287484\n",
       "6     4654.088022   32897.991195  0.011841 -10.248633 -4.436195\n",
       "7     5059.624238   39191.088204  0.008488 -10.437975 -4.769057\n",
       "8     4872.726798   47733.198064  0.009946 -10.665705 -4.610572\n",
       "9     5016.744444   60975.390602  0.004505 -10.932368 -5.402655\n",
       "10    4672.357166  104956.863777  0.005643 -11.515766 -5.177292"
      ]
     },
     "execution_count": 502,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 503,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAELCAYAAAARNxsIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAppUlEQVR4nO3deXhU9dn/8fdN2HfZ1xjWqIgIRqwbiiBYahW1uFuXKvWx1p221Nra2v5qBUTFBSmuT5+2WutuNRAFcUMJgoBi2LeA7DsJ2e7fHzPQBBJIyJkt83ld11yc+c45Z+4Mh3w4y9zH3B0REZEg1Ip1ASIiUnMoVEREJDAKFRERCYxCRUREAqNQERGRwChUREQkMLVjXUCstWrVytPS0mJdhohIQpk9e/Ymd2994HjSh0paWhrZ2dmxLkNEJKGY2cryxnX4S0REAqNQERGRwChUREQkMAoVEREJjEJFREQCk/RXf4mIJJPX5+QyJjOHtdvy6NC8AaOGpjO8b8fA1q9QERFJEq/PyWX0q/PJKywGIHdbHqNfnQ8QWLDo8JeISJIYk5mzP1D2ySssZkxmTmDvoVAREUkSa7flVWn8SChURESSRIfmDao0fiQUKiIiSWLU0HQa1EkpM9agTgqjhqYH9h46US8ikiT2nYzX1V8iIhKI4X07BhoiB9LhLxERCYxCRUREAqNQERGRwCTEORUzOxGYCNQHioBb3P2LcuZbAewEioEid8+IYpkiIkkvIUIFeAj4vbu/a2bDws/PrmDege6+KWqViYjIfoly+MuBpuHpZsDaGNYiIiIVSJQ9lTuATDMbSygIT6tgPgemmJkDT7v7pCjVJyIixFGomFkW0K6cl+4FBgF3uvu/zexS4BlgcDnznu7ua82sDTDVzL519xnlvNdIYCRAampqYD+DiEiyM3ePdQ2HZWbbgebu7mZmwHZ3b3qYZe4Hdrn72EPNl5GR4dnZ2cEVKyKSBMxsdnkXQyXKOZW1wFnh6XOAxQfOYGaNzKzJvmlgCLAgahWKiEj8HP46jJuAR82sNpBP+NCVmXUAJrv7MKAt8FpoR4bawN/d/b0Y1SsikpQSIlTc/WPgpHLG1wLDwtPLgD5RLk1EREpJlMNfIiKSABQqIiISGIWKiIgERqEiIiKBUaiIiEhgFCoiIhIYhYqIiARGoSIiIoFRqIiISGAUKiIiEhiFioiIBEahIiIigVGoiIhIYBQqIiISGIWKiIgERqEiIiKBUaiIiEhgFCoiIkmkpMR5Z946Rkz8lF17iwJff0LcTlhERKrH3Zm+aCNjM3P4eu0OerZtzLptefRo2yTQ91GoiIjUcJ8v28zYKTnMWrGVzi0aMP6yPlzQpyMptSzw91KoiIjUUPPXbGfMlBxmLNpImyb1+OPw47k0ozN1a0fuzIdCRUSkhlm8ficPT13Euwu+o3nDOvx62DH8+NQ06tdJifh7J0SomFkfYCLQGFgBXOXuO8qZ7zzgUSAFmOzuD0azThGRWFq9ZQ+PZC3mtTlraFi3NrcP6sGNZ3ahSf06UashIUIFmAzc4+4fmtkNwCjgvtIzmFkK8ARwLrAGmGVmb7r7N1GvVkQkijbsyOfxaUv4xxerqGXGjWd25eazutGiUd2o15IooZIOzAhPTwUyOSBUgP7AEndfBmBm/wQuBBQqIlIjbd1dwMQZS3nh0xUUFTuXndyZn5/Tg3bN6sespkQJlQXABcAbwAigcznzdARWl3q+Bjgl8qWJiETXrr1FPPvxcv46Yxm7CooYfmJH7hjcg6NbNop1afETKmaWBbQr56V7gRuAx8zst8CbQEF5qyhnzCt4r5HASIDU1NQjqldEJNryC4v528yVPDl9KVt2FzDkuLbcPSSd9HbBftekOuImVNx98GFmGQJgZj2BH5Tz+hrK7sF0AtZW8F6TgEkAGRkZ5QaPiEi8KCwu4ZXZa3js/cWs257PmT1acfeQdE7s3DzWpR0kbkLlUMysjbtvMLNawG8IXQl2oFlADzPrAuQClwNXRrFMEZFAlZQ4b81by8NTF7Fy8x76pjZn3KV9OK1bq1iXVqGECBXgCjP7WXj6VeA5ADPrQOjS4WHuXmRmtxI6iZ8CPOvuX8emXBGRI+fuZC3cwLgpOXz73U6OadeEZ67N4Jxj2mAW/Lfgg2TuyX30JyMjw7Ozs2NdhogIAJ8u2cRDmTnMXb2NtJYNuWtIOuf3bk+tCLRUqQ4zm+3uGQeOJ8qeiohIjTZn1VbGTsnhkyWbad+sPg9e3JtLTupEnZTEaiavUBERiaFvv9vBuCmLmPrNelo2qst95x/HVaekRqWlSiQoVEREYmDFpt2Mz1rEm1+tpXG92twzpCfXn96FRvUS+9dyYlcvIpJg1m3P47H3l/By9mrqpBg3n9WNnw7oSvOG0W+pEgkKFRGRKNi8ay9PTV/KizNX4u5cfUoqPxvYnTZNY9dSJRIUKiIiEbQjv5DJHy3nmY+WkVdYzMX9OnH7oB50btEw1qVFhEJFRCQC8gqKefGzFTz14VK27SlkWO923HVuT7q3iZ+WKpGgUBERCVBBUQkvzVrFhA+WsGHnXs5Ob809Q9I5vmOzWJcWFQoVEakxXp+Ty5jMHNZuy6ND8waMGprO8L4do/Le+YXFHHPfe/ufn5x2FI9f2Y/+XVpE5f3jhUJFRGqE1+fkMvrV+eQVFgOQuy2P0a/OB4hosBSXOL1+9x75hSX7x569LoOB6fHfUiUSFCoiUiOMyczZHyj75BUWMyYzJyKh4u6c98hH5KzfuX+sToqx8A/nUTvBvgUfJIWKiNQIa7flVWm8Om56MZup36wvM7bwD+fRoG5ifgs+SAoVEakROjRvQG45AdKheYPA3uP+N7/m+U9XlBmbc9+5HBWDe8HHK4WKiNQIo4amlzmnAtCgTgqjhqZXe92TZizl//3n2zJjH/9yIJ2OqpnfNakOhYqI1Aj7zpsEefXX63NyueOluWXG3r39TI5t37Q6pdZoChURqTGG9+0YyEn5jxZv5Jpnvigz9vebTonrOy7GC4WKiEjYgtztnD/h4zJjj1/Zl/NP6BCjihKPQkVEkt6qzXsYMGZambHfnn8cN5zRJUYVJS6Fiogkrc279nLSH7PKjI0c0JVfDzs2RhUlPoWKiCSdPQVFHPfbzDJjP+jdnieu6hejimoOhYqIJI3C4hJ63PtumbETOjXjjZ+dnpQtVSJBoSIiNZ67c8ZfppX5cmTT+rWZ+9sh1KqlMAlSQoSKmfUBJgKNgRXAVe6+o5z5VgA7gWKgyN0zolimiMShK/86k0+Xbi4zlvPH86hXWy1VIqHKoWJmjYB8dy8+7MzBmQzc4+4fmtkNwCjgvgrmHejum6JXmojEo1++Mo+XsleXGZt3/xCa1q8To4qSw2FDxcxqAZcDVwEnA3uBema2EfgPMMndF0e0SkgHZoSnpwKZVBwqIpLEHs1azPisRWXGZo4eRLtmNete8PGqMnsq04AsYDSwwN1LAMysBTAQeNDMXnP3v0WuTBYAFwBvACOAzhXM58AUM3PgaXefFMGaRCSOrN6yhwfe/oYppboHZ901oMbfvjfemLsfegazOu5eWN15DluIWRbQrpyX7gVygMeAlsCbwG3u3rKcdXRw97Vm1obQHs3P3X1GOfONBEYCpKamnrRy5crqlC4iMbR+Rz4TPljMS7NWU8uM1k3q8ejlJ3LS0cl1x8VoM7PZ5Z23PmyoHGKFzdx9e7Urq/r79gT+5u79DzPf/cAudx97qPkyMjI8Ozs7wApFJBq27i5g4odLef7TFRSXOJf378zPz+lB26Y6zBUNFYVKda7+esPMNgPfAXOAF929oBrrq5CZtXH3DeHzO78hdCXYgfM0Amq5+87w9BDgD5GoR0RiZ9feIp75aDmTP1rGroIiLjqxI3cM7klqS7WhjwfVCZVP3P3e8C/wMUBP4BfBlHWQK8zsZ+HpV4HnIHS4C5js7sOAtsBr4S8w1Qb+7u7vRageEYmy/MJi/jZzJU9OX8qW3QWc16sddw3pSc+28XHO5PU5uYG23U9U1QmVo8wsA5gHNCX0/ZCIcPdHgUfLGV8LDAtPLwP6RKoGEYmNwuISXs5ezYT3l/DdjnzO7NGKe4ak06dz81iXtt/rc3LL3CAsd1seo1+dD5B0wVKdULkDuCX8eBvQxd8iEpjiEuetr9YyPmsRKzfv4aSjj2L8ZSdyareDrtGJuTGZOWXuOAmQV1jMmMwchUplhc+fPBJcKSIioZYqU79Zz7gpi8hZv5Nj2zfl2esyGJjeJm77c60t1f6lMuM1WUK0aRGR5PDJkk08lJnDV6u30bVVIx6/si/Djm8f9/25OjRvUKavWOnxZBNIqIRbpyxz9+lBrE9EksuXq7YyNjOHT5dupkOz+jx0yQlc3K8jtVNqxbq0Shk1NL3MORWABnVSGDU0PYZVxUZQeypvA30DWpeIJImF63Ywbsoishaup1Xjuvzuh8dxRf9U6tdJrGaP+86b6OqvyvX+Knb3lAPG/gCkAHOBueHeX5nlLC4icpAVm3YzPmsRb361lsb1ajNqaDrXnZZGo3qJe0R+eN+OSRkiB6rM36ABmFljd98F4O6/NbO2hPZOLjGzbu5+UwTrFJEaYO22PCZ8sJiXs9dQN6UW/3NWN346oBvNGuri0ZqiMqGyr4/LHDO7EFjk7kXuvt7MlukLhiJyOJt37eXJ6Uv535krcXeu+d7R3DKwG22aqKVKTVOVfc0WwBPAMWa2CVgI9AO6R6IwEUl8O/ILmTxjGc98vJy8wmJ+dFInbhvUg05HqaVKTVWVUFnl7gMBzKwTcCyQfBdhi8hh5RUU88JnK3hq+lK25xXygxPac+fgnnRv0zjWpUmEVSVUmpvZacDX7r4GWBOhmkQkQRUUlfDPWauY8MESNu7cy8D01tw9JJ3jOzaLdWkSJVUJlUbAPUAvM6tH6PDXAncfFZHKRCRhFJc4r83J5ZGsRazZmkf/tBY8eVU/Tk7TPU2STaWv/gJOdfelAGZWHzgO6BWpwkQk/rk77y34jnFTF7Fkwy56d2zGny7qzYAeraLWUkXdgePLYUPF3WuF/1xaaiwf+DL8EJEk4+7MWLyJsZk5zM/dTvc2jZl4dT+G9moX1f5c6g4cfxL3m0YiEhOzVmxhTGYOXyzfQqejGjBuRB+G9+1ISgz6c6k7cPypcqiY2Q/d/a1IFCMi8WtB7nbGTslhes5GWjepxwMX9uKyk1OpWzt2/bnUHTj+HMmeyp8AhYpIkliyYRfjpy7infnraN6wDr/6/jFce2oaDerGvj+XugPHnyMJlfjuQS0igVizdQ+PZi3m31+uoUGdFG4b1IMbz+xC0/rx01JF3YHjz5GEih9+FhFJVBt25vPktKX83+crMTNuOL0L/3N2N1o2rhfr0g6i7sDxRyfqRQSA7XsKmThjKc9/soKC4hIuzejMbYO6075ZfB9KUnfg+KJQEUlyu/cW8dwny3l6xjJ27S3igj4duHNwT9JaNYp1aZKAjiRU1gdehYhEXX5hMX//fBVPTFvC5t0FDD62LXcP6cmx7ZvGujRJYFUOFXc/NxKFmNkI4H5CjSr7u3t2qddGAz8BioHb3P2gG4KZWQvgJSANWAFc6u5bI1GrSCIrKi7h31+u4dGsxazdns+pXVsy6rx0+qUeFevSpAaIp8NfC4CLgadLD5rZccDlhFrCdACyzKynuxcfsPyvgPfd/UEz+1X4+S8jX7ZIYigpcd6Zv47xUxexbNNuTuzcnDEj+nB691axLk1qkLgJFXdfCJTX4uFC4J/uvhdYbmZLgP7AZ+XMd3Z4+gVgOgoVEdydaTkbGJO5iIXrdpDetgl//XEGg49tE9WWKpIcqhUqZna9uz8XVDEV6AjMLPV8TXjsQG3dfR2Au68zszYVrdDMRgIjAVJTUwMsVSS+fLZ0M2Myv+XLVds4umVDHr38RM4/oQMptUyNGCUiqrun8nug0qFiZllAu3Jeutfd36hosXLGqvVdGXefBEwCyMjI0PdupMb5avU2xk7J4aPFm2jXtD7/76LejMjoRJ2UUEsVNWKUSDlsqJjZvIpeAtpW5c3cfXBV5g9bA3Qu9bwTsLac+dabWfvwXkp7YMMRvJdIQlu0fifjpuSQ+fV6WjSqy29+cCxXf+9o6tcp21JFjRglUiqzp9IWGAoceCWVAZ8GXtHB3gT+bmYPEzpR3wP4ooL5rgUeDP9Z0Z6PSI2zavMeHslaxGtzc2lctzZ3nduTG87oQuN65f8TVyNGiZTKhMrbQGN3n3vgC2Y2PahCzOwiYALQGnjHzOa6+1B3/9rMXga+AYqAn+278svMJgMTw5cfPwi8bGY/AVYBI4KqTSRefbc9nwkfLOalWaupnWKMHNCVmwd046hGdQ+5nBoxSqSYe3KfUsjIyPDs7OzDzygSR7bsLmDih0t54dMVFJc4V/RP5dZzutO2af1KLX/gORUINWL888W9dfhLKsXMZrt7xoHjlTmnYn6Y5KnMPCJSfTvzC3nm4+VM/mg5ewqKGN63I3cM6klqy4ZVWo8aMUqkVObw1zQz+zfwhruv2jdoZnWBMwidv5gGPB+RCkWE/MJi/vezlTw5fQlb9xTy/ePbcde5PenRtskRr1ONGCUSKhMq5wE3AP8wsy7ANqA+kAJMAcaXd75FRKqvsLiEl7NX89j7i1m/Yy8DerbmniE9OaFT81iXJlKuw4aKu+cDTwJPmlkdoBWQ5+7bIlybSNIqLnHe/CqX8VMXs2rLHjKOPopHL+/L97q2jHVpIodUpS8/unshsC5CtYgkPXdnyjfrGTclh0Xrd3Fc+6Y8e10GA9PVUkUSQ6VDxczuKmd4OzBbh79Eqsfd+WRJqKXKV2u207V1I564sh/fP74dtWopTCRxVGVPJSP8eCv8/AfALOBmM/uXuz8UdHEiyWD2yq2Mzczhs2Wb6di8AQ/96AQu7tuR2uGWKiKJpCqh0hLo5+67AMzsd8ArwABgNqBQEamChet2MDYzh/e/3UCrxnW5/4fHccUpqdSrHWqpooaPkoiqEiqpQEGp54XA0e6eZ2Z7gy1LpOZatnEX47MW89ZXa2lavzajhqZz/elpNKz733+OavgoiaoqofJ3YKaZvUGo79f5hC4zbkSohYqIHMLabXk89v5i/jV7DXVTanHrwO7cNKArzRrUOWheNXyURFXpUHH3B8zsP4S+8GjAzaVu+XtVJIoTqQk27drLk9OW8reZKwH48alHc8vZ3WndpF6Fy6jhoySqqt5PpQgoIXQ/k8LgyxGpObbnFTL5o2U88/Fy8guLGXFSZ24b3IOOlWjaqIaPkqiqcknx7cBNwL8J7an8zcwmufuESBUnkoj2FBTx/KcrePrDZWzPK+T8E9pz57k96da6caXXMWpoerkNH0cNTY9EySKBqcqeyk+AU9x9N4CZ/YXQfeIVKiLA3qJi/vnFaiZ8sIRNu/ZyzjFtuHtIT3p1aFbldanhoySqqoSKAaXPHBZT/q1+RZJKUXEJr83J5ZGsxeRuy+OULi14+pp+nHR0i2qtVw0fJRFVJVSeAz43s9cIhclw4NlIFCWSCEpKnPe+/o5xU3JYunE3J3RqxoOX9OaM7q3UUkWSVlWu/no4fKfH0wmFyrVqzyLJyN2Zvmgj46bksCB3Bz3aNGbi1ScxtFdbhYkkvcrcpGsnoau99g+Ves3dvWkkChOJR9c99wXTczYC0LlFAx6+tA8XntiRFPXnEgEq1/r+yO8CJBIlkW5pcu9r8/m/z/ffo46fDujK3UPSqVtb/blESqvq91RE4k4kW5o8MW0JYzJzyox9Nvoc2jfT90VEyqNQkYQXiZYmr8xewz3/+qrM2JQ7B9CzGrfvFUkGcRMqZjYCuB84FuhfqgUMZjaa0PdkioHb3D2znOXvJ/TlzI3hoV+7+38iXLbEgSBbmkzL2cD1z80qM/byT0+lf5fqXR4skiziJlSABcDFwNOlB83sOOByoBfQAcgys57uXnzwKhjv7mMjXqnElSBamny1ehsXPvFJmbGJV/fjvOPbV7s+kWQSN6Hi7guB8i7JvBD4p7vvBZab2RKgP6Fv84tUq6XJ8k27GTh2epmxBy7sxTWnpgVcpUhyiJtQOYSOwMxSz9eEx8pzq5n9GMgG7nb3rZEuTmLvSFqabNy5l5P/lFVm7NaB3blHvbVEqiWqoWJmWUC7cl66193fqGixcsa8nLGngAfCrz0AjANuqKCOkcBIgNTU1MNULYmgsi1Ndu0t4vjflT0ld1Hfjoy/7MQIVSaSXKIaKu4++AgWWwN0LvW8E7C2nHWv3zdtZn8F3j5EHZOASQAZGRnlBZTUMAVFJfT8zbtlxjKOPopX/ue0GFUkUjMlwuGvN4G/m9nDhE7U9wC+OHAmM2vv7uvCTy8idOJfklxJiXPKn99n487/3vG6TZN6zBw9iFr6FrxI4OImVMzsIkJt9FsD75jZXHcf6u5fm9nLhG5ZXAT8bN+VX2Y2GZgYvvz4ITM7kdDhrxXAT2PwY0gcGTHxU2atKHtabdEfv69vwYtEkLkn99GfjIwMz87OPvyMkjDuemkur87JLTO24PdDaVwvbv4PJZLwzGy2u2ccOK5/ZVJjzFy2mcsnzSwz9sW9g2jTpH6MKhJJPgoVSXgrN+9m/NRFvD73v9dvTLvnbLq0ahTDqkSSk0JFEtZ32/N57IPFvDxrNbVTjJvP6sbNZ3WlecO6EXvPSHdDFkl0ChVJOFt2F/DU9CW8+NlKSty58pRUbh3YnTZNI3uYK5LdkEVqCoWKJIyd+YVM/mg5kz9aRl5hMRf368Ttg3rQuUXDqLx/JLohi9Q0ChWJe3kFxbz42Qqe+nAp2/YUMqx3O+46tyfd20S3DX2Q3ZBFaiqFisStgqISXspezYT3F7Nh517O6tmae4ak07tTs5jUE0Q3ZJGaTqEicae4xHljbi7jsxaxekseJ6cdxYQr+nJK15Yxras63ZBFkoVCReKGu5P59XrGTclh8YZd9OrQlOeuP56ze7Yu75YIUXck3ZBFko1CRWLO3flo8SbGTslh3prtdG3diCeu7Mf3j28Xd/25KtsNWSRZKVQkpmav3MJD7+Xw+fItdGzegDE/OoGL+nakdor6c4kkIoWKxMTXa7czbsoiPvh2A60a1+P3F/Ti8v6dqVc7JdaliUg1KFQkqpZt3MXDUxfx9rx1NK1fm1+cl851p6XRsK42RZGaQP+SJSpyt+XxWNZiXvlyDfVq1+LWgd25aUBXmjWoE+vSRCRAChWJqI079/Lk9CX838xVAFx7ahq3DOxGq8b1YlyZiESCQkUiYnteIX+dsYxnP1nO3qISRpzUiZ8P6kFHfVFQpEZTqEig9hQU8dwnK3j6w6XsyC/ih306cOfgHnRt3TjWpYlIFChUJBB7i4r5x+ereHzaUjbt2sugY9pw15Ce9OoQm5YqIhIbChWplqLiEl6dk8ujWYvJ3ZbHKV1a8PQ1/Tjp6BaxLk1EYkChcgR0oyYoKXHeXfAd46bmsGzjbk7o1IwHL+nNGd1bxUVLFRGJDYVKFSX7jZrcnek5GxmTmcM363bQs21jJl59EkN7tVWYiIhCpaqS+UZNny/bzJjMHLJXbqVziwaMv6wPF/TpSEqc9ecSkdiJm1AxsxHA/cCxQH93zw6PtwReAU4Gnnf3WytYvgXwEpAGrAAudfetQdeZjDdqmr9mO2Om5DBj0UbaNKnHH4cfz6UZnalbW/25RKSsuAkVYAFwMfD0AeP5wH3A8eFHRX4FvO/uD5rZr8LPfxl0kcl0o6bF63fy8NRFvLvgO5o3rMOvhx3Dj09No34d9ecSkfLFTai4+0LgoOPy7r4b+NjMuh9mFRcCZ4enXwCmE4FQSYYbNa3esodHshbz2pw1NKxbm9sH9eDGM7vQpL5aqojIocVNqASgrbuvA3D3dWbWJhJvUpNv1LRhRz4TPljCP2etopYZN57ZlZvP6kaLRnVjXZqIJIiohoqZZQHtynnpXnd/I4p1jARGAqSmplZ5+Zp2o6atuwuYOGMpL3y6gqJi57KTO/Pzc3rQrln9WJcmIgkmqqHi7oMjuPr1ZtY+vJfSHthwiDomAZMAMjIyPII1xbVde4t49uPl/HXGMnYVFDH8xI7cMbgHR7dsFOvSRCRB1aTDX28C1wIPhv+M2p5PoskvLOZvM1fy5PSlbNldwJDj2nL3kHTS2zWJdWkikuDiJlTM7CJgAtAaeMfM5rr70PBrK4CmQF0zGw4McfdvzGwyMDF8+fGDwMtm9hNgFTAiBj9GXCssLuFf2Wt47P3FfLcjnzO6t+Keoemc2Ll5rEsTkRoibkLF3V8DXqvgtbQKxm8sNb0ZGBSR4hJcSYnz1ry1PDx1ESs376FvanMevqwPp3VrFevSRKSGiZtQkeC5O1kLNzBuSg7ffreTY9o14ZlrMzjnmDZqqSIiEaFQqaE+XbKJhzJzmLt6G2ktG/LYFX05v3d7aqmliohEkEKlhpmzaitjp+TwyZLNtG9Wnwcv7s0lJ3WiTopaqohI5ClUaog35uZy50tzKXFo2agu951/HFedkqqWKiISVQqVBDc9ZwPXPTdr//O7zu3JT87oQqN6+qsVkejTb54ENWfVVi568tMyY9PvOZu0VvrioojEjkIlwSzZsJPBD88oM/b2z8/g+I66F7yIxJ5CJUGs3ZbHaQ9+UGbsHzd9j1O7tYxRRSIiB1OoxLktuwsY/PCHbNldsH9s4tX9OO/49jGsSkSkfAqVOLV7bxEXPfkJi9bv2j/254t7c0X/qndVFhGJFoVKnCkoKuH657/gkyWb94+NGprOzwYe7h5lIiKxp1CJEyUlzl0vz+X1uWv3j11/ehq/Pf84tVQRkYShUIkxd+dP7yxk8sfL94+df0J7Hr28LylqqSIiCUahEkMTP1zKg+9+u//5KV1a8OJP+lOvtr4FLyKJSaESAy9nr+YXr8zb/7xb60a8cesZNNa34EUkwem3WBRN/WY9N72Yvf95swZ1+ODus2jZuF4MqxIRCY5CJQq+WL6FS5/+rMzYx78cSKejGsaoIhGRyFCoRNDCdTv4/qMflRnLvGOA7gUvIjWWQiUCVm3ew4Ax08qMvXLzqWSktYhRRSIi0aFQCdDGnXs5a8w09hQU7x977rqTGXhMmxhWJSISPQqVAOzML+T8CR+zcvOe/WMPX9qHi/t1imFVIiLRp1CphvzCYq6e/DnZK7fuH/vND47lxjO7xrAqEZHYiZsbl5vZCDP72sxKzCyj1HhLM5tmZrvM7PFDLH+/meWa2dzwY1gk6/3P/HUcc997+wPl5rO6sfzPwxQoIpLU4mlPZQFwMfD0AeP5wH3A8eHHoYx397ERqO0gi8Pdgy/p14kxPzqBWmqpIiISP6Hi7guBg5onuvtu4GMzi6s2vbcP7sHtg3vEugwAXp+Ty5jMHNZuy6ND8waMGprO8L4dY12WiCShuDn8FZBbzWyemT1rZkfFuphoeH1OLqNfnU/utjwcyN2Wx+hX5/P6nNxYlyYiSSiqoWJmWWa2oJzHhQGs/imgG3AisA4Yd4g6RppZtpllb9y4MYC3jp0xmTnkFRaXGcsrLGZMZk6MKhKRZBbVw1/uPjiC616/b9rM/gq8fYh5JwGTADIyMjxSNUXD2m15VRoXEYmkGnP4y8xK37T9IkIn/mu8Ds0bVGlcRCSS4iZUzOwiM1sDnAq8Y2aZpV5bATwMXGdma8zsuPD45FKXHz9kZvPNbB4wELgzuj9BbIwamk6DOmXvv9KgTgqjhqbHqCIRSWbxdPXXa8BrFbyWVsH4jaWmr4lMZfFt31VeuvpLROJB3ISKHLnhfTsqREQkLsTN4S8REUl8ChUREQmMQkVERAKjUBERkcCYe0J/96/azGwjsDLWdVSgFbAp1kUcguqrHtVXPaqveqpb39Hu3vrAwaQPlXhmZtnunnH4OWND9VWP6qse1Vc9kapPh79ERCQwChUREQmMQiW+TYp1AYeh+qpH9VWP6queiNSncyoiIhIY7amIiEhgFCoiIhIYhUqMmdn9ZpZrZnPDj2HlzNPZzKaZ2UIz+9rMbq/K8pGuLzzfeWaWY2ZLzOxXpcZbmNlUM1sc/jMit3k2s3vMzM2sVTmvpZeqf66Z7TCzO6ry80WyvvDrK8K3bphrZtmlxuPh84vZ9leZ+sKvx2T7M7MHwrcwn2tmU8ysQznzxGz7q0x94fmC2/7cXY8YPoD7gXsOM097oF94ugmwCDiusstHob4UYCnQFagLfFWqvoeAX4WnfwX8JQI1dgYyCX2JtVUlav2O0Be3Iv75VbY+YEV5r8XD5xfL7a+S9cVs+wOalpq+DZgYT9tfZesLcvvTnkoCcPd17v5leHonsBCIp173/YEl7r7M3QuAfwIXhl+7EHghPP0CMDwC7z8e+AVQmatOBgFL3T2aXRSqUt+BYv75xcH2d7jPL2bbn7vvKPW00SFq3Ceq298R1HegKn9+CpX4cGt4F/XZw+1emlka0Bf4/EiWj1B9HYHVpZ6v4b+/dNq6+zoI/XIC2gRZmJldAOS6+1eVXORy4B8HjEXs86tCfQ5MMbPZZjay1HhcfX7R3v4qWV/Mtr9wjX8ys9XAVcBvDzN7VLe/KtQX3PYXyd1WPfbvQmYBC8p5XAi0JbRLXAv4E/DsIdbTGJgNXFxqrNLLR6o+YAQwudTza4AJ4eltB8y7NeD6PgeahedbwSEOfxE6NLIp/A8lWp9fpeoDOoT/bEPo8M2AOPz8YrH9Hba+WG5/B8w3Gvh9PG1/la0vyO2vSsXrEdkHkAYsqOC1OoSOK991JMtHsj7gVCCz1PPRwOjwdA7QPjzdHsgJsJ7ewIbwL5sVQBGwCmhXwfwXAlOi9flVtb5Sy91P+Dh7vHx+sdj+KltfrLa/cuo9+lA/f7S3v6rWF9T2p8NfMWZm7Us9vYjQ/zAOnMeAZ4CF7v5wVZePdH3ALKCHmXUxs7qEdvHfDL/2JnBtePpa4I2ganP3+e7ext3T3D2N0GGPfu7+XQWLXMEBhx4i+flVtj4za2RmTfZNA0NK1RHzzy9W218V/n5jsv0BmFmPUk8vAL49xOxR3f7C6z9sfYFvf5FKRT0q/b+H/wXmA/PCf4H7/lfQAfhPePoMQsc85wFzw49hh1o+mvWFnw8jdFXQUuDeUuMtgfeBxeE/W0Tws1xB+PBIOfU1BDYTPpRyuJ8vmvURumrpq/Dj63j7/GK5/VXh7zcm2x/wb0K/gOcBbwEd42n7q0x9QW9/atMiIiKB0eEvEREJjEJFREQCo1AREZHAKFRERCQwChUREQmMQkVERAKjUBE5DDPbVY1lG5jZh2aWEmRNsWJmdc1shpnVjnUtEp8UKiKRdQPwqrsXx7qQIHioC/D7wGWxrkXik0JFpArM7C4zWxB+3FFq/D4z+zZ8I6N/mNk94ZeuItzawsx6m9knpZbpZ2YfRLjef5nZ42b2sZmtNLMzzOxFM1tkZs8c4WpfJ/RziRxEu7AilWRmJwHXA6cABnxuZh8S6jJ7CaGW8LWBL4HZ4T5UXd19RXgVXwPdzCwlvOcyDrg7wmX3Bj5z91vN7A+EenidTahlyHozu8Xd91ZxnQuAk4MtU2oKhYokPTPLAtqV89K97l66gd4ZwGvuvju83KvAmYT2+N9w97zw+Fvh+VsB2/Yt7O4lZvY10Cvc6G+Vh29+ZWYPuPt9QdZsZvWB5sAj4dfygGc8fH8MM9sDFFT1/d292MwKzKyJh27aJbKfQkWSnrsPruSsVsXxPKD+AWMzgdOBW4DzAMysHQf8WzSzeoTCYGt4/uHuvrWKNfcCvnT3kvDzPsBT4fV3Ata6u5f3/pWooR6QX4kaJMnonIpI5c0AhptZw3CL8IuAj4CPgR+aWX0zawz8ACD8CzglvMewz0zgj4T2eHLDY30Jdf4t7RbgOXf/NbCldKBUQW9CnWf3OYFQt1oIBcy+6fLev8IazKwlsNHdC4+gJqnhFCoilRQ+VPU88AWhOxJOdvc57j6LUNvyr4BXgWxge3ixKYQOm+3zLbAX+EupsRM5+Jd6H2BeOKQquj/M4fTet95wsDUoFU6lA6a89z9UDQOB/xxhTVLDqfW9SADMrLG77zKzhoT2aEa6+5dm1pfQ3RKvCc/3ODDL3V8otewzwE1Aa+B8d3/GzK4Azgd2AHPcfVIEaz/o/cPj5dYQPpc02t1zIlWTJC7tqYgEY5KZzSV05de/952Ad/c5wDQz62Zm3xLaW3ih9ILu/pPweY++wPLwcB1gLbAbeDGShVfw/uXWEL6i7XUFilREeyoiIhIY7amIiEhgFCoiIhIYhYqIiARGoSIiIoFRqIiISGAUKiIiEhiFioiIBEahIiIigVGoiIhIYP4/FjHH3a2paUQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.scatter(table['rhs'],table['lhs'])\n",
    "m, b = np.polyfit(table['rhs'], table['lhs'], 1)\n",
    "plt.plot(table['rhs'],table['rhs']*m + b)\n",
    "plt.xlabel('$-\\log(y_{i,g}-m_{i,g})$')\n",
    "plt.ylabel('$\\log(1-\\overline{\\pi}_{i,g})$')\n",
    "plt.savefig('../figures/fig_a1_production.eps')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 504,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/users/loulou/anaconda3/lib/python3.8/site-packages/scipy/stats/_stats_py.py:1477: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=10\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>           <td>lhs</td>       <th>  R-squared:         </th> <td>   0.846</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.826</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   43.87</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 27 May 2022</td> <th>  Prob (F-statistic):</th> <td>0.000165</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:38:12</td>     <th>  Log-Likelihood:    </th> <td> -3.3449</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    10</td>      <th>  AIC:               </th> <td>   10.69</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>     8</td>      <th>  BIC:               </th> <td>   11.29</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -4.1207</td> <td>    0.908</td> <td>   -4.537</td> <td> 0.002</td> <td>   -6.215</td> <td>   -2.026</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>rhs</th>       <td>    1.3589</td> <td>    0.205</td> <td>    6.624</td> <td> 0.000</td> <td>    0.886</td> <td>    1.832</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.662</td> <th>  Durbin-Watson:     </th> <td>   1.956</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.436</td> <th>  Jarque-Bera (JB):  </th> <td>   1.050</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.537</td> <th>  Prob(JB):          </th> <td>   0.591</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 1.831</td> <th>  Cond. No.          </th> <td>    35.3</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                    lhs   R-squared:                       0.846\n",
       "Model:                            OLS   Adj. R-squared:                  0.826\n",
       "Method:                 Least Squares   F-statistic:                     43.87\n",
       "Date:                Fri, 27 May 2022   Prob (F-statistic):           0.000165\n",
       "Time:                        15:38:12   Log-Likelihood:                -3.3449\n",
       "No. Observations:                  10   AIC:                             10.69\n",
       "Df Residuals:                       8   BIC:                             11.29\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -4.1207      0.908     -4.537      0.002      -6.215      -2.026\n",
       "rhs            1.3589      0.205      6.624      0.000       0.886       1.832\n",
       "==============================================================================\n",
       "Omnibus:                        1.662   Durbin-Watson:                   1.956\n",
       "Prob(Omnibus):                  0.436   Jarque-Bera (JB):                1.050\n",
       "Skew:                           0.537   Prob(JB):                        0.591\n",
       "Kurtosis:                       1.831   Cond. No.                         35.3\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 504,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = ols('lhs ~ rhs',data=table).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 495,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.822957452418385"
      ]
     },
     "execution_count": 495,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(model.params[1]-0.5)/model.HC0_se[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 496,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.36222965765550713"
      ]
     },
     "execution_count": 496,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/model.params[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 497,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07759032081961857"
      ]
     },
     "execution_count": 497,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(1/(model.params[1]**2))*model.HC0_se[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 498,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/users/loulou/anaconda3/lib/python3.8/site-packages/scipy/stats/_stats_py.py:1477: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=10\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>           <td>lhs</td>       <th>  R-squared:         </th> <td>   0.728</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.650</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   9.375</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 27 May 2022</td> <th>  Prob (F-statistic):</th>  <td>0.0105</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:37:35</td>     <th>  Log-Likelihood:    </th> <td> -7.0821</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    10</td>      <th>  AIC:               </th> <td>   20.16</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>     7</td>      <th>  BIC:               </th> <td>   21.07</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   15.8955</td> <td>   21.412</td> <td>    0.742</td> <td> 0.482</td> <td>  -34.736</td> <td>   66.527</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>rhs</th>       <td>   11.5976</td> <td>   11.915</td> <td>    0.973</td> <td> 0.363</td> <td>  -16.576</td> <td>   39.771</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>rhs2</th>      <td>    1.2271</td> <td>    1.652</td> <td>    0.743</td> <td> 0.482</td> <td>   -2.679</td> <td>    5.133</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 0.589</td> <th>  Durbin-Watson:     </th> <td>   1.247</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.745</td> <th>  Jarque-Bera (JB):  </th> <td>   0.555</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.226</td> <th>  Prob(JB):          </th> <td>   0.758</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 1.939</td> <th>  Cond. No.          </th> <td>1.86e+03</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.86e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                    lhs   R-squared:                       0.728\n",
       "Model:                            OLS   Adj. R-squared:                  0.650\n",
       "Method:                 Least Squares   F-statistic:                     9.375\n",
       "Date:                Fri, 27 May 2022   Prob (F-statistic):             0.0105\n",
       "Time:                        15:37:35   Log-Likelihood:                -7.0821\n",
       "No. Observations:                  10   AIC:                             20.16\n",
       "Df Residuals:                       7   BIC:                             21.07\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     15.8955     21.412      0.742      0.482     -34.736      66.527\n",
       "rhs           11.5976     11.915      0.973      0.363     -16.576      39.771\n",
       "rhs2           1.2271      1.652      0.743      0.482      -2.679       5.133\n",
       "==============================================================================\n",
       "Omnibus:                        0.589   Durbin-Watson:                   1.247\n",
       "Prob(Omnibus):                  0.745   Jarque-Bera (JB):                0.555\n",
       "Skew:                          -0.226   Prob(JB):                        0.758\n",
       "Kurtosis:                       1.939   Cond. No.                     1.86e+03\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.86e+03. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 498,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table['rhs2'] = table['rhs']**2\n",
    "model = ols('lhs ~ rhs + rhs2',data=table).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "c06ef71ca6a83554d8b678a29f56d641660e46430c130240f7802f9a0facbb6c"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
